© 2025. The Author(s). Published by the American Astronomical Society.
Nonlinear Coupling between Magnetar QPOs
P. Pnigouras1	and S.K. Lander2
Departamento de Física Aplicada, Universidad de Alicante, Campus de San Vicente del Raspeig, Alicante E-03690, Spain
2 School of Engineering, Mathematics and Physics, University of East Anglia, Norwich NR4 7TJ, UK
Received 2025 January 21; revised 2025 April 7; accepted 2025 April 17; published 2025 June 25
Abstract
The quasiperiodic oscillations (QPOs) observed in the tails of magnetar giant γ-ray flares have long been interpreted as normal oscillation modes of these stars. However, most studies modeling QPOs have neglected some key features in the analyses of the signals, namely that QPOs appear to be detectable only intermittently and exhibit drifts in their frequencies. These are typical characteristics of nonlinear mode coupling, where, at leading order, the modes couple and evolve collectively as triplets. Using a representative triplet of modes, we solve the system’s nonlinear equations of motion analytically and argue that the coupling is likely axial-axial-polar in nature, with the observed intermittence and frequency drifts providing a way to infer details of the magnetar’s internal magnetic field geometry.
Unified Astronomy Thesaurus concepts: Magnetars (992); Stellar oscillations (1617); Soft gamma-ray repeaters (1471); Pulsation modes (1309); Stellar pulsations (1625); Neutron stars (1108); Magnetic fields (994); Stellar magnetic fields (1610); Neutron star cores (1107); Analytical mathematics (38); Polytropes (1281); Non-radial pulsations (1117)

1.		Introduction
Two decades ago, quasiperiodic oscillations (QPOs) were discovered (G. L. Israel et al. 2005; T. E. Strohmayer & A. L. Watts 2005) in the aftermath of giant γ-ray flares from two different magnetars, neutron stars (NSs) with very strong magnetic fields B ∼ 1014–1015 G characterized by their high- energy outburst activity. These QPOs, in the frequency range of tens to thousands of Hz, seem to represent the first detections of oscillation modes of isolated NSs and so should encode valuable information about their interior physics, including the critically important NS equation of state (A. W. Steiner & A. L. Watts 2009), for which there is substantial uncertainty at high densities.
Although the lower-frequency QPOs had the expected properties of elastic shear modes of the NS crust, it quickly became clear that the problem was far more complex. The crust and core are threaded by a strong magnetic field that couples their dynamics, but whose broad structure is not well understood; this introduces an unwanted degree of ambiguity in the interpretation of observed QPOs as certain modes of the star. Another particularly concerning issue for QPO interpretation has been the possibility that the core admits a continuum of local oscillation frequencies rather than global oscillation modes, as this continuum would tend to drain energy from typical crustal modes (Y. Levin 2007; M. van Hoven & Y. Levin 2012). Although certain additional pieces of physics appear to remove or alter the continuum—e.g., axisymmetric modes with a radial component (H. Sotani & K. D. Kokkotas 2009; A. Colaiuda & K. D. Kokkotas 2012), nonaxisymmetric modes (S. K. Lander et al. 2010; S. K. Lander & D. I. Jones 2011), or the inclusion of neutron superfluidity (M. Gabler et al. 2013b; A. Passamonti & S. K. Lander 2013, 2014)—these also further complicate
Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. attempts to identify particular observed QPOs with theoretically predicted modes.
Alongside modeling challenges, there have also been various reanalyses of the original data, resulting in conflicting conclusions about how significant and long-lived the different QPOs are (V. Hambaryan et al. 2011; D. Pumpe et al. 2018; M. C. Miller et al. 2019). Although a few additional individual QPOs have been detected from various magnetars (e.g., A. J. Castro-Tirado et al. 2021, O. J. Roberts et al. 2023), these are no substitute for one new set of high-quality data that reveals multiple QPOs, and modeling has reached an impasse of simply waiting passively for the next giant flare.
Almost all work on modeling QPOs has focused on matching numbers: tuning the parameters of a magnetar model so that its oscillation-mode frequencies agree with observed QPOs. Unfortunately, the theoretically predicted modes are sufficiently densely spaced in frequency that comparison with observations does not place robust constraints on the magnetar’s underlying physics. Furthermore, this body of work ignores a major feature of the observations: the QPOs do not behave as if they were all excited at once and then simply decay, like the ringing of a bell. Instead, over hundreds of seconds following the giant flare, QPOs appear and disappear, suggesting some process of decay and reexcitation.
In this paper, we propose that these features are due to nonlinear couplings between different modes. Modeling the couplings would not only help understand the full oscillation tails rather than just single modes, but could also lead to new insight about the physical origin of these modes.
2.		Revisiting the QPOs
To study magnetar mode coupling, we need to identify potential candidate QPOs that could be exhibiting this behavior. We will look for QPOs following the giant flares of soft γ-ray repeaters (SGR) 1806-20 and 1900+14 (hereafter SGR 1806 and SGR 1900) that were detected for more than just a very short segment, and preferably that are common to both magnetars and have been corroborated by more than one study. For SGR 1806, these criteria lead us to select the QPOs at 18, 26, 30, 93, and 150 Hz; there are also a few intriguing QPOs in the kHz range (G. L. Israel et al. 2005; T. E. Strohmayer & A. L. Watts 2006). A subsequent search for QPOs from the far more frequent smaller bursts by D. Huppenkothen et al. (2014) revealed an additional QPO at 57 Hz, adding weight to the 59 Hz QPO found by V. Hambaryan et al. (2011) in an earlier analysis of the giant-flare tail from the same magnetar. For SGR 1900, the original data are of lower quality, but the analysis of T. E. Strohmayer & A. L. Watts (2005) nonetheless yielded QPOs at 28, 53, 84, and 155 Hz.
Since mode coupling is known to cause frequency drift (see Section 5), we assume that QPOs of similar frequency represent the same underlying mode; i.e., for SGR 1806, we regard the strong 92.5 Hz and weak 95 Hz QPOs of G. L. Israel et al. (2005) and the 90 Hz QPO of T. E. Strohmayer & A. L. Watts (2006) as the same mode; likewise for the 30 Hz (G. L. Israel et al. 2005) and 26Hz (A. L. Watts & T. E. Strohmayer 2006) QPOs. This leaves us with four modes observed from both magnetars, with frequencies ∼28, 55, 88, and 150 Hz. We will disregard the ∼150 Hz mode, partly because it is a very broad peak and partly because it is closer to the frequency range at which the core continuum of frequencies may affect results.
Our ansatz is that each of the above QPOs represents an oscillation mode of the star of magnetoelastic character. To identify candidate physical modes, we use the fitting formulae from M. Gabler et al. (2016, Equation (20)). For the modes to have such a character and not be damped into the continuum, it is required that 1014 G  B  5 × 1015 G. We can also exclude any oscillations with nonconstant phase, whose eigenfunction does not reach the stellar surface, from our modeling.
Satisfying the above, the only plausible triplet involves three of the “lUn modes” of M. Gabler et al. (2016), namely
2U2, 6U2, 6U4. (1)
Varying the field strength so that the frequencies best match observations, we infer B = 8.6 × 1014 G for SGR 1806, which gives triplet mode frequencies of
29.8,		56.2, 92.3 Hz, (2)
and B = 7.2 × 1014 G for SGR 1900, yielding
27.3,		53.7, 84.0 Hz. (3)
These values of B are in very good agreement with those inferred from spindown measurements: 7 × 1014 G for SGR 1900 (S. Mereghetti et al. 2006) and 8 × 1014 G for SGR 1806 (C. Kouveliotou et al. 1998; G. Younes et al. 2017). Note that the latter is often reported to have a substantially stronger field, 2 × 1015 G, but this value was inferred during a period of enhanced spindown around the time of the giant flare (D. M. Palmer et al. 2005)—likely due to an extended corona, making the dipole-spindown formula unreliable.
3.		Nonlinear Mode-coupling Theory
Nonlinear mode coupling is a result of the nonlinear interaction among the oscillation modes of a system. Ever since the work of W. Dziembowski (1982), this mechanism has been considerably studied in the context of astrophysics, mainly to explain the phenomenology of some variable stars (e.g., Y. Wu & P. Goldreich 2001), to study the saturation of NS instabilities (A. K. Schenk et al. 2001; P. Arras et al. 2003; P. Pnigouras & K. D. Kokkotas 2015, 2016), and to model the effects of nonlinear tides in binary systems (N. N. Weinberg et al. 2012, 2013).
The most basic system is that of three modes, which couple to each other at second order (see, e.g., A. H. Nayfeh & D. T. Mook 1979). In the absence of dissipation, the equations of motion for this system may be written in terms of the mode amplitudes Qα:
Q a= iwakQbQge-iDwt,	(4a)
Qb= iwbkQg*QaeiDwt,	(4b)
Qg= iwgkQaQb*eiDwt,	(4c)
where ωα is the mode frequency, Δω = ωα - ωβ - ωγ is the detuning parameter, and κ is the mode-coupling coefficient, defined as an overlap integral among the eigenfunctions xia of the three modes (for an explicit expression of κ see, e.g., A. K. Schenk et al. 2001 or W. Dziembowski 1982). The mode amplitudes Qα are defined by the decomposition of a generic perturbation ξ i in terms of the modes as
xi(xi, t) = å Qa(t)xia(xi)eiwat	(5)
a
(note that Qα does not contain the mode phase).
In spherical, i.e., nonrotating and unmagnetized, stars, each mode can be identified by a single spherical harmonic Ylm (q, f) in a vector spherical harmonic decomposition, namely
xi =WYlm ir+ViYlm -iUijkjYlmkr.	(6) r
Each mode is characterized by a multipole degree l,an azimuthal order m, and its overtone n, denoting the number of nodes in the mode’s radial profile, described by the functions W (r), V(r), and U(r). According to the standard mode classification, we have polar modes, in which U = 0, and axial modes, in which W = V = 0.
One may then show that
kµ ∬ Ya*Yb Yg sin qdqdf,	(7) with Ya º Ylma. This is a known integral (e.g., J. J. Sakurai & J. Napolitano 2017), which implies that κ ≠ 0if
ma=mb+mg	(8)
and
la=lb+lc-2l,	(9)

where la …lb …lc and l =0, 1,¼lmax  lc/2, with the indices a, b, c taking the values α, β, γ, so that the mode a has the largest degree and the mode c has the lowest. Equations (8) and (9) constitute the selection rules that the coupled mode triplet has to satisfy and restrict the possible couplings. Furthermore, the equations of motion (Equations (4a–4c)) imply that, among the many possible couplings, the triplets that ultimately affect the system’s dynamics the most are those that satisfy the nonlinear resonance condition
Dw=wa-wb-wg 0. (10)
Note that the mode-coupling formalism laid out above is a result of second-order perturbation theory, meaning that these effects could, in principle, manifest themselves in any nonlinear numerical simulation. Apart from offering insight into the otherwise complex nonlinear dynamics, its main advantage is that it does not need to keep track of the mode phase, but rather of the slowly varying (due to nonlinear effects) mode amplitude Q, thus being able to achieve much longer integration times. This is necessary in any attempt to model the QPO features observed over long timescales.
Even though the potential relevance of nonlinear effects for QPOs has been suggested in the past (T. E. Strohmayer & A. L. Watts 2006), previous efforts employing quadratic perturbations were not motivated by the same phenomena that interest us here. For instance, H. Sotani et al. (2024) derived the quadratically perturbed relativistic Euler equation of axisymmetric torsional (axial) modes in the NS crust and subsequently solved it numerically in the time domain, but their main result was to show the indirect excitation of normal modes, as well as the existence of additional modes, in the oscillation spectrum at second order, rather than attempting an explanation of the observed intermittence and frequency drifts of magnetar QPOs. Nonetheless, their calculation constitutes an important step toward a more realistic treatment of the problem and is worth building upon.
4.		Modeling the Best Triplet

xf(r, q) = Csin a

(2n+ 1)p r
2R*

¶qYl0, sinq

(11)

The magnetoelastic modes of M. Gabler et al. (2016) depend on the magnetic field structure within the star, with nodes and maxima following the magnetic field lines. Thus, their indices n and l (see Equation (1)) do not have the same meaning as above. An accurate description of such a mode would therefore require an infinite sum of Ylm components in a vector spherical harmonic decomposition. This would present severe difficulties in a modecoupling calculation, due to the various multipolar contributions from each mode entering the coupling coefficient (Equation (7)). However, aiming for quantitatively accurate mode amplitudes in a calculation that is otherwise necessarily qualitative would be somewhat pointless. Instead, we argue that it is the similarity of the mode eigenfunctions involved in the nonlinear coupling that matters most, effectively establishing how good the overlap among them is, and hence, how strongly they are coupled. With this in mind, we will assume that the coupling of three magnetic lUn modes is comparable to the coupling of three modes with eigenfunctions of the form of Equation (6) with the same l, n, and therefore proceed by representing each relevant mode from M. Gabler et al. (2016) as where R* is the stellar radius and C is a constant determined by the chosen normalization. In the above we have adopted the form of an axisymmetric (m = 0) axial mode, in accordance with M. Gabler et al. (2016), and have also assumed that the mode’s radial dependence is sinusoidal (see, e.g., U. Lee 2008), with the requisite number of nodes n, and nonzero at the surface.
Like almost all other calculations of magnetar QPOs (e.g., L. Samuelsson & N. Andersson 2007; H. Sotani et al. 2007;
P. Cerdá-Durán et al. 2009), the work of M. Gabler et al. (2016) assumes axial oscillations, since these are more readily excited: they do not involve a radial displacement, which is inhibited by the composition-gradient stratification of the star (R. C. Duncan 1998). We therefore begin with the result of A. K. Schenk et al. (2001), who calculated the coupling coefficient for a triplet consisting of two axial modes β, γ, and a third mode α whose Eulerian density perturbation δρ vanishes:
kabg =12òp(G1-G)kxak* i xbj jxigdV,(12) with G = dlnp/dlnr and G1 =(¶lnp/¶lnr)xp, where p is the pressure and ρ is the density. The adiabatic exponent Γ1 differs from the background index Γ if the chemical composition (here denoted by the proton fraction xp) changes throughout the star.
Since the divergence of an axial mode is zero, Equation (12) shows that there is no coupling between a triplet of axial modes. Thus, mode α must have a polar component. It may still have a leading-order axial piece, however; for instance, as is known from stellar perturbation theory, rotation introduces polar corrections to axial modes, and vice versa (e.g., W. Unno et al. 1989). Therefore, only the polar piece of mode α contributes to ixia. We will regard this divergence as a parameter to adjust; from Equation (6) we have ixia =fa Ya, where
fa =rz2[¶r(rWa)-la(la+1)Va], (13) with the small parameter ζ quantifying the relative amplitude of the polar piece to the axial piece of the mode.
There are various possible mechanisms for axial-polar mixing, but perhaps the most compelling is related to the magnetic field geometry. The purely poloidal background magnetic fields employed in M. Gabler et al. (2016) and most other magnetar mode modeling are dynamically unstable (P. Markey & R. J. Tayler 1973; G. A. E. Wright 1973), meaning that only a linked poloidal–toroidal geometry is realistic. Such a field geometry leads to axial modes mixing with polar modes in a manner proportional to the relative strength of toroidal to poloidal field (A. Colaiuda & K. D. Kokkotas 2012). We will therefore regard ζ as a measure of the strength of this toroidal component.
Using Equation (6) and assuming axisymmetry, the axial- axial-polar coupling coefficient (12) takes the form
kabg =4∬ Ya*Yb Yg sin qdqdf ò p (G1 -G)
´ ⎧⎨ fa(rgb¶rgg+r¶rgbgg)(-La+Lb+Lg)
⎡ (La + Lb + Lg)2 222 + fagbgg⎢ 2 -(La+Lb+Lg)
+ 3(-La + Lb + Lg)⎥⎬ r2dr,
(14) where Λ = l(l + 1) and g = U/r2. For the evaluation of the functions W, V, and U, we will assume that the polar components ξr and ξθ of mode α take on the same radial
dependence as the axial component ξ f, according to the ansatz (11).
We also define the mode energy as Ea= 2Ia w2a∣Qa∣2, where

I = ò xi*xirdV

(15)

(see, e.g., A. K. Schenk et al. 2001). The modes will be normalized by fixing their energy to unit amplitude, i.e., 2Ia wa2 = Mc2, with M being the mass of the star and c being the speed of light. Hence, the mode amplitude is simply related to the mode energy through Eα = |Qα|2Mc2.
The same normalization has been tacitly assumed for the coupling coefficient κ in the derivation of Equations (4a–4c); κ is calculated in units of the chosen energy value—in this case Mc2—so the expression of Equation (14) needs to be divided by that before being used in Equations (4a–4c). An additional subtlety is that the modes β and γ also ought to have polar components of order ζ, which would thus couple to the axial components of the other two modes in the triplet (we will, however, neglect couplings of order ζ2 and higher, involving the polar components of two or three modes). This means that we also need to evaluate kb¯ga¯ and kg¯a¯b, with the bar over the mode label meaning that the mode eigenfunction has to be complex conjugated in Equation (14). The complete expression for the coupling coefficient is then
1
k=	(kabg +kb¯ ga¯+kg¯a¯b).
Mc2

(16)

We will assume a background star with M = 1.4 Me (Me being the solar mass) and R* = 10 km, described by a polytropic equation of state p ∝ ρΓ, with Γ = 2, for which the background quantities may be obtained analytically, and set the adiabatic exponent Γ1 = 2.1 (see, for example, N. Andersson & P. Pnigouras 2020).
5.		Results
The coupled equations of motion (Equations (4a–4c)) admit an analytical solution. Decomposing the complex mode amplitude Qα into its real amplitude and phase components as
iJa
Qa = eae	,	(17) k wbwg
Equations (4a–4c) become
e a= ebegsin j,	(18a)
eb=-egeasin j,	(18b)
eg=-eaebsinj,	(18c) and
j =cosj⎜⎛ebeg-egea-eaeb⎟⎞ +Dw,	(18d) ⎝ ea	eb	eg ⎠
where j = ϑα - ϑβ - ϑγ + Δωt.
Including the phase ϑα in the harmonic time dependence of the mode, we get xia µ ei(wat+Ja), suggesting that the nonlinear coupling induces modulation of the mode frequency of the ¢
form wa =wa+Ja, with the shift Ja given by
ebeg
Ja =	cos j.	(19) ea

of integration:	
1 = e2a + e2b,	(20)	
22 2 =ea+eg,	(21)	
and		
Dw 2 L =eaebegcos j+ea , 2a	(22)	

Upon inspection, from the equations of motion

(Equations (18a–18d)) we can obtain the following constants

with 1 and 2 being related to the energy of the mode triplet, whereas L dictates the evolution of the frequency shifts Ja. Combining Equations (18a) and (22), making the replacement ea2 = q1, we get

1	Dw	L 2⎤
q(1-q)(m-q)-1⎜⎛⎝	2 q-	1⎟⎞⎠	⎦⎥,

⎛dq⎞2=4 ⎝ dt ⎠	1 (23) with m = 2/1. The expression in brackets on the right-hand side of Equation (23) is a cubic polynomial with roots q1, q2, and q3, such that q ä [q1, q2 „ 1] and q3 …1.
We further replace q = q1 + (q2 - q1)sin2 d, which makes Equation (23)
⎛ dddt ⎞ = 1(q3 - q1)(1 - k2 sin2 d),	(24) with k2 = (q2 - q1)/(q3 - q1) „ 1. We integrate to get
d	dd¢ ò	=1(q3-q1)t,(25) d0	1- k2 sin2 d¢ from which we can define the Jacobi elliptic function sn(u, k)= sin d (e.g., M. Abramowitz & I. A. Stegun 1972), where d0	dd
u(t)=1(q3-q1)t+ò	,	(26) 0	1- k2sin2d with the angle δ0 determined by the initial conditions for the mode amplitudes.
Switching to our initial variables, the solution to Equations (4a–4c) is
2 ∣Qa∣	=1	[q1+(q2-q1)sn2(u, k)],	(27a)	
	wb wg k2		
∣Qb∣2 =	1
1	[1 -q1-(q2-q1)sn2(u, k)],	(27b)	
	wg wa k2		
2 ∣Qg∣ =	1	[m - q1 - (q2 - q1)sn2(u, k)].	(27c)	
	wa wb k2		

These are periodic with respect to time, with the period given by
2K (k)
1(q3 - q1)
K(k) being the complete elliptic integral of the first kind.
To demonstrate the system’s behavior, we will use the mode triplet of SGR 1900 (Equation (3)), with α labeling the high- frequency (parent) mode and β, γ the mid- and low-frequency

Figure 1. Evolution of our representative coupled mode triplet, with α labeling the high-frequency (parent) mode and β, γ the mid- and low-frequency (daughter) modes, respectively. We plot the mode amplitude |Q| (top) and mode frequency shift Dnº J/2p in Hz (bottom) as functions of time, for ζ = 4 × 10-4 and initial conditions |Qα| = 10-3, |Qβ| = |Qγ| = 10-5.

mag a

L q2- q1
1 q1q2 ,

Dw
2

mag Jb

 mag
g

(29a)

(29b)

(29c)

(daughter) modes, respectively. The results are similar for the mode triplet in SGR 1806 (Equation (2)). The initial mode amplitudes |Q| are fixed, leaving ζ as the only free parameter.
In Figure 1 we plot the evolution of the mode amplitudes |Q| and frequency shifts Dnº J /2p (in Hz) for ζ = 4 × 10-4 and initial conditions |Qα| = 10-3, |Qβ| = |Qγ| = 10-5. The energy exchange between the modes, characteristic for these systems, occurs periodically (T ≈ 0.7 s), with the high-energy mode pumping energy into the low-energy modes, which then transfer it back, and the cycle repeats. Likewise, the frequency shift also oscillates periodically, with its maxima corresponding to minima in the amplitude and reaching values up to 4 Hz.
The period T of the amplitude modulation and the magnitude of the frequency shifts can be adjusted by changing ζ. Unless ζ falls below a certain limit value zlim (see below), (q3 - q1) ∼ 1 and K(k) ∼ 1, so from Equation (28) we have T ~ 1-1/2 µ k-1 µ z-1. The scaling of Δν with ζ is more complicated. Expressing Equation (22) as L =ea2 (Ja + Dw/2) and noticing that ea2 Ja =e2bJb =e2gJg, one may derive the magnitudes of the frequency shifts as
L	q2- q1
1 (1-q1)(1-q2),
Dw L	q2- q1
2	1 (m-q1)(m-q2)
These are increasing functions of ζ and approach zero in the linear limit, when	(or q ≈ q ; see below). Hence, ,	lim 1	2	.	,
Pnigouras & Lander decreasing ζ leads both to larger modulation periods and smaller frequency shifts.
There is, however, a limit value, zlim, to which we alluded earlier, below which the modes are too weakly coupled and their amplitudes are practically constant (essentially, the system approaches the linear limit). This occurs when q1 ≈ q2, or when the discriminant of the cubic polynomial in q from Equation (23) gets close to zero. The value of zlim depends on the value of the detuning Δω: the larger Δω is, the larger zlim is. This is sensible, because a small detuning implies a stronger nonlinear resonance and, hence, stronger nonlinear coupling among the three modes (see Equation (10)), which allows for smaller values of the coupling coefficient.
For the triplet we are examining, the case depicted in Figure 1 is already near zlim . In order to explain the appearance and disappearance of modes over timescales of order 10–100 s, which is what QPO observations suggest, we need lower values of ζ. Let us, therefore, reduce the original value of the detuning (Δω/2π = 3Hz) by a factor of 1000. For this case, we plot |Q| and Δν in Figure 2 for three different values of ζ using the same initial conditions as before, in order to demonstrate that the variation of ζ affects T and Δν according to the scalings presented before. In particular, for ζ = 5 × 10-6 we obtain T ≈ 50 s, and Δνα reaches values up to 1 Hz; for ζ = 10-5 the period drops to T ≈ 25 s, and Δνα increases to up to 3 Hz; and for ζ = 3 × 10-5 the period is T ≈ 9 s, with Δνα getting up to 30 Hz.
Moreover, in Figure 3 we plot |Q| and Δν for the same parameters as in the first column of Figure 2, but also assuming an amplitude detection cutoff, set to |Q| = 4 × 10-4. When the amplitude becomes lower than this cutoff value, the mode is regarded as undetectable, and its amplitude and frequency shift curves are plotted as semitransparent. This way we mimic the observed appearance and disappearance of modes in giant-flare tails.
Finally, note that our choice of initial conditions is such that zlim  1 . Increasing the initial mode amplitudes practically has the same effect as increasing ζ (and vice versa), which is to be expected; a large mode amplitude would imply that we are deep into the nonlinear regime, hence the coupling coefficient (parameterized by ζ) need not be as large for the characteristic energy exchange to occur among the three modes.
6.		Discussion
The discovery of magnetar QPOs triggered a concerted effort to develop theoretical models of oscillating NSs, with the intention of identifying the observed QPOs as specific oscillation modes of the model magnetar. This work was, ultimately, motivated by the prospect of using QPO identification as a way to probe unknown aspects of the interior physics of magnetars: their magnetic field strength, crustal properties, and the equation of state of matter beyond nuclear density (A. W. Steiner & A. L. Watts 2009; H. Sotani et al. 2012; M. Gabler et al. 2013a). Although this endeavor has not been as fruitful as first hoped, due both to the limited quality of the data and to the complexity of the theoretical modeling problem, the advances made to date will be invaluable in the interpretation of QPOs following any future giant flare. Furthermore, this work has since found a separate application, related to QPOs associated with γ-ray bursts (GRBs). One class of GRBs is believed to be related to compact-object binary mergers, and QPOs were recently detected in a GRB precursor event

Figure 2. The effect of varying ζ on |Q| (top panels) and Δν (bottom panels), for the same initial conditions as for Figure 1, but with the detuning dropped by a factor of 1000. From left to right: ζ = 5 × 10-6, 10-5, 3 × 10-5.

Figure 3. Evolution of |Q| and Δν from the first column of Figure 2, but now assuming an amplitude detection cutoff |Q| = 4 × 10-4 (dashed line).

(C. Chirenti et al. 2024; S. Xiao et al. 2024), providing insights into the premerger physics (A. G. Suvorov et al. 2022).
The two sets of magnetar-QPO data represent more than just a set of numbers to be matched, however. In this paper we revisit the original X-ray tails following the giant flares of SGR 1806 and SGR 1900, arguing that these hold additional valuable information that deserves study. In particular, the apparent disappearance and reappearance of magnetar QPOs, together with their drifting frequencies, are important observations that have previously received little attention from theoretical modeling (a rare exception being the work of Y. Levin 2007). We have found that both phenomena can be explained through nonlinear mode coupling. In such a system, mode amplitudes have oscillatory variation in time (rather than simply decaying), and also time-varying frequency shifts, providing a concrete physical reason to believe that, e.g., the 26 and 30 Hz QPOs of SGR 1806 are the same oscillation mode. Despite earlier recognition of the potential importance of nonlinear effects in QPO modeling, e.g., by H. Sotani et al. (2024), where the second-order axial oscillation spectrum of the NS crust was derived, their capability of explaining these QPO features had not been appreciated before.
There is one key free parameter ζ in our coupled mode triplet, which quantifies the strength of the subdominant polar piece of a predominantly axial mode. This is required, because the most readily excited magnetar modes are axial, but three purely axial modes do not experience any coupling. The value of ζ dictates the timescale on which a given mode dis-/ reappears, and also its drift in frequency. We argue that it is natural to associate ζ with the toroidal component of the star’s magnetic field; the behavior of a magnetar’s QPOs could thus be one of very few possible mechanisms for inferring the star’s internal field geometry.
Although we can mimic QPO dis-/reappearance over the observed timescales, or observed shifts in frequency of some QPOs, it is more challenging to get both right with our simple coupled-triplet model—partly because we predict the maximum frequency shift to occur when the mode is least likely to be detectable. Furthermore, the real data do not indicate that QPO variability is strictly periodic, as our model predicts. It is, however, easy to see where further complexity could come from in the mode coupling, perhaps resulting in better agreement with observations. First, the fact that the ∼57, ∼90, and ∼150 Hz QPOs are also likely to couple together (based on the frequency selection rule, i.e., 57 + 90 ≈ 150)—involving two of the same modes as our fiducial triplet here—suggests a possible coupling between the two triplets. Further complexity would come if we allowed for energy to be dissipated from the modes (rather than just transferred between them). Finally, we need to recall that we do not observe these QPOs directly, as oscillations of the star’s surface, but rather as modulations to the brightness of a fireball anchored to the star by magnetic field lines (M. Gabler et al. 2014; T. Akgün et al. 2018). This process will affect how “visible” certain modes are, on top of the intrinsic changes in amplitude that mode coupling gives.
Higher-frequency QPOs (600 Hz) are more likely to represent nonaxisymmetric modes, whose eigenfunctions would be sufficiently different as to inhibit coupling with axisymmetric mode triplets like those we focus on here. In the tail following the giant flare of SGR 1806, these high-frequency QPOs exhibited similar behavior to the lower-frequency QPOs (A. L. Watts & T. E. Strohmayer 2006), suggesting a process of mode coupling in this case too. Intriguingly, high-frequency (kHz) QPOs were also seen during the giant flare from an extragalactic magnetar (A. J. Castro-Tirado et al. 2021), providing an additional probe of the underlying physics driving this phenomenon (H. Sotani 2024) and suggesting that these kHz oscillations may be preferentially excited during the initial phase, perhaps coupling to other modes much later on.
Acknowledgments
P.P. acknowledges support from the María Zambrano Fellowship Programme (ZAMBRANO21-22), funded by the Spanish Ministry of Universities and the University of Alicante through the European Union’s “Next Generation EU” package, as well as from the grant PID2021-127495NB-I00, funded by MCIN/AEI/ 10.13039/501100011033 and by the European Union, from the Astrophysics and High Energy Physics program of the General- itat Valenciana ASFAE/2022/026, funded by the Spanish Ministry of Science and Innovation (MCIN) and the European Union’s “Next Generation EU” package (PRTR-C17.I1), and from the Prometeo 2023 excellence program grant CIPROM/ 2022/13, funded by the Ministry of Education, Culture, Universities, and Occupation of the Generalitat Valenciana.
ORCID iDs
P. Pnigouras	https: /orcid.org/0000-0003-1895-9431
References
Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover),
Akgün, T., Cerdá-Durán, P., Miralles, J. A., & Pons, J. A. 2018, MNRAS, 481, 5331
Andersson, N., & Pnigouras, P. 2020, PhRvD, 101, 083001
Arras, P., Flanagan, É. É., Morsink, S. M., et al. 2003, ApJ, 591, 1129 Castro-Tirado, A. J., Østgaard, N., Göğüş, E., et al. 2021, Natur, 600, 621
Cerdá-Durán, P., Stergioulas, N., & Font, J. A. 2009, MNRAS, 397, 1607
Chirenti, C., Dichiara, S., Lien, A., & Miller, M. C. 2024, ApJ, 967, 26
Colaiuda, A., & Kokkotas, K. D. 2012, MNRAS, 423, 811
Duncan, R. C. 1998, ApJL, 498, L45
Dziembowski, W. 1982, AcA, 32, 147
Gabler, M., Cerdá-Durán, P., Font, J. A., Müller, E., & Stergioulas, N. 2013a, MNRAS, 430, 1811
Gabler, M., Cerdá-Durán, P., Stergioulas, N., Font, J. A., & Müller, E. 2013b, PhRvL, 111, 211102
Gabler, M., Cerdá-Durán, P., Stergioulas, N., Font, J. A., & Müller, E. 2014, MNRAS, 443, 1416
Gabler, M., Cerdá-Durán, P., Stergioulas, N., Font, J. A., & Müller, E. 2016, MNRAS, 460, 4242
Hambaryan, V., Neuhäuser, R., & Kokkotas, K. D. 2011, A&A, 528, A45
Huppenkothen, D., Heil, L. M., Watts, A. L., & Göğüş, E. 2014, ApJ, 795, 114
Israel, G. L., Belloni, T., Stella, L., et al. 2005, ApJL, 628, L53
Kouveliotou, C., Dieters, S., Strohmayer, T., et al. 1998, Natur, 393, 235
Lander, S. K., & Jones, D. I. 2011, MNRAS, 412, 1730
Lander, S. K., Jones, D. I., & Passamonti, A. 2010, MNRAS, 405, 318
Lee, U. 2008, MNRAS, 385, 2069
Levin, Y. 2007, MNRAS, 377, 159
Markey, P., & Tayler, R. J. 1973, MNRAS, 163, 77
Mereghetti, S., Esposito, P., Tiengo, A., et al. 2006, ApJ, 653, 1423
Miller, M. C., Chirenti, C., & Strohmayer, T. E. 2019, ApJ, 871, 95
Nayfeh, A. H., & Mook, D. T. 1979, Nonlinear Oscillations (New York: John Wiley & Sons)
Palmer, D. M., Barthelmy, S., Gehrels, N., et al. 2005, Natur, 434, 1107
Passamonti, A., & Lander, S. K. 2013, MNRAS, 429, 767
Passamonti, A., & Lander, S. K. 2014, MNRAS, 438, 156
Pnigouras, P., & Kokkotas, K. D. 2015, PhRvD, 92, 084018
Pnigouras, P., & Kokkotas, K. D. 2016, PhRvD, 94, 024053
Pumpe, D., Gabler, M., Steininger, T., & Enßlin, T. A. 2018, A&A, 610, A61
Roberts, O. J., Baring, M. G., Huppenkothen, D., et al. 2023, ApJL, 956, L27
Sakurai, J. J., & Napolitano, J. 2017, Modern Quantum Mechanics (2nd ed.;
Cambridge: Cambridge Univ. Press)
Samuelsson, L., & Andersson, N. 2007, MNRAS, 374, 256
Schenk, A. K., Arras, P., Flanagan, É. É., Teukolsky, S. A., & Wasserman, I.
2001, PhRvD, 65, 024001
Sotani, H. 2024, PhRvD, 109, 023030
Sotani, H., & Kokkotas, K. D. 2009, MNRAS, 395, 1163
Sotani, H., Kokkotas, K. D., & Stergioulas, N. 2007, MNRAS, 375, 261
Sotani, H., Nakazato, K., Iida, K., & Oyamatsu, K. 2012, PhRvL, 108, 201101
Sotani, H., Suvorov, A. G., & Kokkotas, K. D. 2024, PhRvD, 110, 023035
Steiner, A. W., & Watts, A. L. 2009, PhRvL, 103, 181101
Strohmayer, T. E., & Watts, A. L. 2005, ApJ, 632, L111
Strohmayer, T. E., & Watts, A. L. 2006, ApJ, 653, 593
Suvorov, A. G., Kuan, H. J., & Kokkotas, K. D. 2022, A&A, 664, A177
Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial
Oscillations of Stars (2nd ed.; Tokyo: Univ. Tokyo Press)
van Hoven, M., & Levin, Y. 2012, MNRAS, 420, 3035
Watts, A. L., & Strohmayer, T. E. 2006, ApJL, 637, L117
Weinberg, N. N., Arras, P., & Burkart, J. 2013, ApJ, 769, 121
Weinberg, N. N., Arras, P., Quataert, E., & Burkart, J. 2012, ApJ, 751, 136
Wright, G. A. E. 1973, MNRAS, 162, 339
Wu, Y., & Goldreich, P. 2001, ApJ, 546, 469
Xiao, S., Zhang, Y.-Q., Zhu, Z.-P., et al. 2024, ApJ, 970, 6
Younes, G., Baring, M. G., Kouveliotou, C., et al. 2017, ApJ, 851, 17

M. Gabler et al. (2016) classification (see Section 2).