Couplings of torsional and shear oscillations in a neutron star crust

Hajime Sotani [email protected] Astrophysical Big Bang Laboratory, RIKEN, Saitama 351-0198, Japan Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), RIKEN, Saitama 351-0198, Japan Theoretical Astrophysics, IAAT, University of Tübingen, 72076 Tübingen, Germany Arthur G. Suvorov Theoretical Astrophysics, IAAT, University of Tübingen, 72076 Tübingen, Germany Departament de Física Aplicada, Universitat d’Alacant, Ap. Correus 99, E-03080 Alacant, Spain Manly Astrophysics, 15/41-42 East Esplanade, Manly, NSW 2095, Australia Kostas D. Kokkotas Theoretical Astrophysics, IAAT, University of Tübingen, 72076 Tübingen, Germany

Abstract

Mature neutron stars are thought to be sufficiently cold that nuclei in the outer layers freeze, solidifying a crust. Crustal elasticity allows the star to support a set of seismic modes, such as torsional oscillations. These axial-parity modes can couple to the polar sector in a number of ways, for example via rotation or a magnetic field. Even in a static, spherically-symmetric star however these modes can couple at non-linear order. In this study, such couplings in the crust are examined for the first time: we derive the axisymmetric perturbation equations for second-order axial eigenfunctions, which are sourced by axial-polar couplings at first-order, and solve the resulting equations in the time domain. Through our studies, we find that the second-order spectrum contains additional oscillation modes, not predicted by the linear analysis with either axial or polar perturbations, which can be excited to relatively large amplitudes.

pacs:

04.40.Dg, 97.10.Sj, 04.30.-w

I Introduction

Neutron stars, produced in supernova explosions ending the lives of massive stars, are some of the most suitable objects to probe physics at extreme scales. Their cores readily reach densities exceeding that of nuclear saturation while maintaining baryonic degrees of freedom at low temperatures, placing them in a parameter space inaccessible to terrestrial laboratories [1]. The elastic, solid crust of a neutron star is also thought to host various exotic states of matter, such as hadron-quark mixed “pasta” phases [2, 3] and ion lattices soaked in superfluid neutrons [4, 5]. With the exception of gravitational waves, the crust effectively mediates all of the persistent and transient activity from these stars, and thus detailed modelling of both dynamical processes and the crust itself is necessary to fully utilize existing and upcoming instruments.

The crust can sustain a variety of seismic oscillation modes [6, 7, 8, 9]. Overstraining events, particularly in the magnetar class of objects where an ultrastrong field is capable of exerting great stress [10, 11, 12], might excite crust-localized [13, 14] and global magnetoelastic [15, 16] modes. This theoretical prediction is supported by observations of quasi-periodic oscillations (QPOs) with frequencies on the order of tens of Hz and up in the X-ray tails of magnetar flares, as in the 1979 event in SGR 0526-66 [17], the 1998 event in SGR 1900+14 [18], and the 2004 event in SGR 1806-20 [19, 20]. The latter demonstrated an especially rich spectrum; studies coming out over a decade later are still finding new Fourier peaks [21, 22, 23]. In principle, the inverse problem of deducing the stellar properties — magnetic geometry, equation of state (EOS), crustal shear modulus, … — can be tackled by matching the observed spectra to detailed models. Such an approach has been carried out by various authors [25, 26, 24, 27, 28, 29, 30], who matched the forest of low- and high-frequency QPOs from various flares to the modes of a general-relativistic star, with a realistic EOS and topologically-mixed, globally-threading magnetic field, exhibiting coupled axial and polar perturbations¹¹1 We note that the axial and polar perturbations are classified by their parity. That is, considering a coordinate transform of polar angles $\theta\to\pi-\theta$ and $\phi\to\pi+\phi$ , an axial quantity of angular order $\ell$ transforms as $(-1)^{\ell+1}$ while polar ones change as $(-1)^{\ell}$ ..

Despite these successes, there remains a number of theoretical challenges of varying severity. Arguably the most serious of these relates to how the failure mechanism actually operates in the crust. The seed for the main event is the evolution of the magnetic field, which builds up mechanical stresses that culminate in a failure event and the release of magnetoelastic energy [31, 32]. This localized relaxation generates horizontal currents that may go on to launch chains of Hall and/or thermoplastic waves that instigate further failures of various amplitude away from the initial site, twisting the external field and priming it for explosive reconnection [33, 34, 35, 36]. Alternatively, the rapid release of poloidal energy that fuels the flare could trigger instability and magnetic reconfiguration [12], exerting Maxwell stress over short timescales. This may help to explain why the dynamic spectra of flare QPOs are so complicated: different regions of the crust fail at different times, leading to the rising and falling of various mode families over the course of hundreds of seconds. This means that a fully-fledged attempt at a solution of the aforementioned inverse problem requires not only the specification of initial data at ignition, but also injections of energy at later times, the particulars of which require detailed, microphysical calculations of failure propagation that are presently out of reach.

A second, related problem is that most studies of magnetar oscillations work within the linear regime (though see [37, 38, 39, 40, 41, 42] for some exceptions). That is, one solves for a small perturbation over a fixed (or dynamical) background in either the time or frequency domain. While this provides a reasonable, leading picture for the spectrum of possible excitations, it does not allow one to study the physics related to saturation via parametric and other instabilities, which are important since the observed QPOs have significances that vary by orders of magnitude [20, 21, 23]; the respective Bayes factors are presumably tied to the amplitude that individual modes achieve. Aside from predicting saturation amplitudes, nonlinear models are needed to deduce which daughter modes can be excited from some initial combination of parents. An example of this is the so-called $p$ - $g$ instability, where nonresonant couplings between these polar modes could instigate orbital energy losses in a binary inspiral, possibly leading to significant dephasing in a gravitational waveform [43]. In the case of QPOs in flare tails, there could be instances of frequency drifting also. For the SGR 1806-20 giant flare, Miller et al. [23] further report a $\approx 28.1$ Hz oscillation starting at $\approx 193$ s post flare and subsequently disappearing, with a new peak at $\approx 27.9$ Hz emerging some $\gtrsim 20$ s later. Field evolution via Hall avalanches or backreaction from excited modes [44], nonlinear effects both, could theoretically cause such a drift (see also Ref. [16]).

In this paper, which we intend to be the first in a series, we make some additional steps toward the nonlinear modelling of crustal oscillations in a neutron star. We derive equations at (Eulerian) second-order for axial perturbations in general relativity with the Cowling approximation, including new expressions for the shear and strain tensors, building on the linear formalism of Schumaker and Thorne [6]. The second-order perturbations are sourced by first-order, coupled polar-axial perturbations, and the equations are solved for a realistic EOS in the time domain. Although nonlinear torsional oscillations have been studied previously (notably in Ref. [41]), sourcing by axial-polar primaries has not. Our scheme is further capable of handling seed injections at later times, to simulate delayed and non-local failings within the crust, and is stable over long (many seconds) timescales. While we take a step back from the realistic picture by not considering crust-core couplings or magnetic fields, assuming that the stellar magnetic field is not so strong that one can neglect such couplings, i.e., the field strength is $\ll 10^{15}$ G [45, 46], certain key effects, which are supported by analytic estimates, are observed in the simulations. For instance, we predict the existence of new pairs of daughter modes that arise from axial-polar couplings at linear order, the amplitude of which is set by that of the parents. Higher quantum-number modes and overtones from the linear spectrum are also excited with a complicated pattern. The main goal here is to derive the relevant equations and solve them in a simplified setting to pave the way for more realistic investigations in forthcoming work.

This paper is organized as follows. In Sec. II, we briefly introduce the (tabular EOS) equilibrium neutron star model adopted in this study and the shear modulus characterizing the crust elasticity. In Sec. III, we derive the two-dimensional perturbation equations for 2nd order axial oscillations through an order-matching procedure, where the 2nd order perturbation equation has a source term composed of the product of linear axial and polar perturbations. The boundary conditions imposed at the both edges of the (elastic) crust region together with the initial conditions are also described. Then, we show several examples of mode excitation by solving, in the time domain, the derived equation in Sec. IV. Finally, in Sec. V we offer some conclusions. In this study, we adopt geometric units with $c=G=1$ for the speed of light, $c$ , and the gravitational constant, $G$ . The metric signature is $(-,+,+,+)$ .

II Equilibrium models

Since we are ultimately interested in developing a scheme that applies to Galactic magnetars, which spin very slowly, the metric and fluid variables can be taken as static to a good approximation. The background neutron star model is further taken to be spherically symmetric and strain-free, whose line element is given in Schwarzschild-like coordinates $\{t,r,\theta,\phi\}$ through

ds^{2}=-e^{2\Phi}dt^{2}+e^{2\Lambda}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta d\phi^{2}\right).

(1)

With this expression, the only non-zero component of the four-velocity is in the $t$ -direction, and is given by

u^{t}=e^{-\Phi}.

(2)

The metric variables $\Phi$ and $\Lambda$ appearing in Eq. (1) together with the fluid variables, $p$ (pressure) and $\varepsilon$ (energy density), are constructed by integrating the Tolman-Oppenheimer-Volkoff equations, assuming an approximate but appropriate EOS for neutron star matter, i.e., zero-temperature matter in beta-equilibrium. In this study, we specifically adopt the SLy4, which is a Skyrme-type EOS with a self-consistent crustal region [47], and focus on stars with mass $M=1.41M_{\odot}$ and radius $R=11.68$ km.

Since torsional and shear oscillations can be excited only inside regions with elasticity, their frequencies depend on the density at the boundaries. The outer boundary, which is the interface between the crust and the envelope (or “ocean”), depends strongly on the physical temperature, meaning a zero-temperature EOS cannot apply [48]. In this study, we simply assume that the transition (baryon) density at this interface is $10^{10}$ g/cm³ and the density matching to the stellar atmosphere is $10^{6}$ g/cm³. Meanwhile, the inner crust-core boundary is also complicated. In a realistic neutron star, it is thought that the crust region is composed not only of phases of the more-familiar spherical nuclei but also of non-spherical nuclei – the so-called pasta phases. Nevertheless, the thickness of the pasta phase is rather narrow, compared to that of spherical nuclei, e.g., [49]. So, in this study, we only consider the phase of spherical nuclei, by assuming that the transition (baryon) density between the crust and core is $1.045\times 10^{14}$ g/cm³. With these transition densities, the radii of the core and the annulus separating the crust, $R_{c}$ , and the envelope, $R_{e}$ , become $R_{c}=10.85$ km and $R_{e}=11.57$ km, respectively.

The shear modulus, $\mu$ , is another important quantity associated with elastically-supported oscillations. Even though there really are multiple shear moduli which depend on the shape of nuclei [50], in this study we neglect the presence of pasta and assume a single contributor. For this spherical-nuclei shear modulus we simply adopt the standard, low-temperature expression proposed in Ref. [51], i.e.,

\mu=0.1194\frac{n_{i}(Ze)^{2}}{a},

(3)

where $e$ , $n_{i}$ , $Z$ , and $a$ respectively denote the elementary charge, ion number density, charge number, and Wigner-Seitz cell radius (i.e., $4\pi a^{3}/3=1/n_{i}$ ).

II.1 Remarks on nonlinear elasticity

The elastic solid that is the neutron star crust is in reality a complicated, chemically- and magnetically-stratified medium with non-trivial composition and temperature gradients [7]. At linear order though, the crust (or any other solid) necessarily abides by a Hookean relationship; that is, the material is such that strain is linearly proportional to stress. The elastic stress-tensor is itself proportional to the Lagrangian displacements at linear order in a linear mode problem. At nonlinear order, however, the Hookean approximation becomes an additional assumption of the model – some kind of “elastic EOS” is needed to close the system.

As described in Section 12 of Ref. [52] (see also [53]), more general elasticity relationships can be explored. A set of tensors can be introduced on the (Lagrangian) matter space, from which a set of scalar invariants can be constructed which represent a type of basis of possible relationships between the shear and strain tensors. These individual terms are weighted by general-relativistic generalizations of the Lamé coefficients. Although an amount of literature is devoted to the study of the shear and bulk modulii, higher-order coefficients have not been examined in a neutron star crust. By including these additional terms though, a more general stress tensor could be built (Eq. (12.19) in Ref. [52]), which will adjust the spectrum.

To make contact with linear studies, however, we adopt a Hookean approximation throughout. Exploration of the non-linear spectra for different stress-strain relationships will be considered elsewhere.

III Perturbation equations

In this study, we adopt the Cowling approximation, i.e., the metric perturbations are neglected so that the matter fields oscillate on a fixed spacetime. This is a reasonable approximation since even the most energetic of polar-parity modes excited in a magnetar flare are unlikely to produce strong gravitational-wave signals [37, 39, 40]. The torsional oscillations, which are classified as axial perturbations and form our primary interest, are less energetic by orders of magnitude.

We further restrict ourselves to axisymmetric perturbations. This approximation is less physically justifiable, especially since axisymmetric states seem to be an evolutionarily-unlikely outcome in magnetars [54, 55], but it is sufficient to demonstrate the main effects of polar-axial coupling at non-linear order. In general, the spatial components²²2Throughout this work, purely spatial components are written with Latin scripts $(i,j,k,\ldots)$ , while spacetime indices are denoted by members of the Greek alphabet $(\alpha,\beta,\gamma,\ldots)$ . of the Lagrangian displacement vector can thus be given by

\xi^{i}=\left(rW(t,r,\theta),V(t,r,\theta),\frac{{\cal Y}(t,r,\theta)}{\sin\theta}\right),

(4)

where $W$ and $V$ are of polar parity and ${\cal Y}$ encompasses axial perturbations. In what follows, perturbed variables associated with such a Lagrangian displacement are written with a prepended $\delta$ , which in general encompasses a hierarchy of terms. For example, the perturbed velocity can be expanded as

\delta u_{\mu}=\delta u_{\mu}^{(1)}+\delta u_{\mu}^{(2)}+\cdots,\\

(5)

where the superscript, $(i)$ , denotes the $i$ -th order perturbations of each quantity³³3We remark that at this point the scheme differs from several previous studies looking at saturation amplitudes [64, 59, 60], not only because of the relativistic framework but because we explicitly include a second-order velocity perturbation, which is often set to zero..

Using the displacement given by Eq. (4), the spatial parts of the perturbed four-velocity can be written as

	$\displaystyle\delta u^{i}$	$\displaystyle=$	$\displaystyle\left(re^{-\Phi}\partial_{t}W,e^{-\Phi}\partial_{t}V,e^{-\Phi}\frac{\partial_{t}{\cal Y}}{\sin\theta}\right),$		(6)
	$\displaystyle\delta u_{i}$	$\displaystyle=$	$\displaystyle\left(re^{-\Phi+2\Lambda}\partial_{t}W,r^{2}e^{-\Phi}\partial_{t}V,r^{2}e^{-\Phi}\partial_{t}{\cal Y}\sin\theta\right),$		(7)

from which one can show that

\displaystyle\delta u^{\alpha}_{\ ;\alpha}

\displaystyle=re^{-\Phi}\partial_{tr}W+re^{-\Phi}\left(\Lambda_{,r}+\frac{3}{r}\right)\partial_{t}W+e^{-\Phi}\partial_{t\theta}V+\frac{\cos\theta}{\sin\theta}e^{-\Phi}\partial_{t}V.

(8)

We note that $\delta u^{\alpha}_{\ ;\alpha}$ is composed of only polar perturbations and $\delta u^{\alpha}_{\ ;\alpha}=0$ even for non-axisymmetric axial perturbations.

Due to the spherically-symmetric background, linear perturbations with axial parity (torsional $t$ - oscillations) are completely decoupled from those with polar parity (interface $i$ -, shear $s$ -, fundamental $f$ -, and pressure $p$ -modes). However, if one considers higher-order perturbations, the axial and polar sectors inevitably couple via non-linear effects. In this study, as a first step, we examine the simplest case with non-linear effects, i.e., 2nd-order perturbations for the torsional oscillations induced by the coupling between the linear perturbations of both axial and polar parity. We note that one can similarly consider the 2nd-order polar perturbations, which are also induced by linear polar and axial perturbations. But, in this study, we try to understand the magnetar QPOs via crustal oscillations (or at least subspectra due to such components), where most effort has gone into axial modes (again cf. crust-core couplings and global modes [15]). Meanwhile, the problem is a bit more involved for polar modes, so we leave it for future work (along with magnetic fields etc), and we focus on 2nd-order axial perturbations in this study.

The total energy-momentum tensor, $T_{\mu\nu}$ , is composed of two parts: a fluid part, $T_{\mu\nu}^{\rm(F)}$ , and a shear part, $T_{\mu\nu}^{\rm(S)}$ . One can derive the fluid part for the $\phi$ -component of the perturbed energy-momentum conservation law, which contributes to the torsional oscillations, as

	$\displaystyle\delta T^{{\rm(F)}\phi\nu}_{\ \ \ \ \ \ ;\nu}=$	$\displaystyle\frac{(\varepsilon+p)e^{-2\Phi}}{\sin\theta}\bigg{\{}\left[1+\frac{\delta\varepsilon+\delta p}{\varepsilon+p}\right]\partial_{tt}{\cal Y}+r(\partial_{t}W)(\partial_{tr}{\cal Y})+(\partial_{t}V)(\partial_{t\theta}{\cal Y})+\left(r\partial_{tr}W+\partial_{t\theta}V\right)(\partial_{t}{\cal Y})$
		$\displaystyle+\frac{\partial_{t}(\delta\varepsilon+\delta p)}{\varepsilon+p}\partial_{t}{\cal Y}+\left[\frac{\varepsilon^{\prime}+p^{\prime}}{\varepsilon+p}-\Phi_{,r}+\Lambda_{,r}+\frac{5}{r}\right]r(\partial_{t}W)(\partial_{t}{\cal Y})+\frac{2\cos\theta}{\sin\theta}(\partial_{t}V)(\partial_{t}{\cal Y})\bigg{\}}\,,$		(9)

where we omit third- and higher-order terms, such as $\delta\varepsilon(\partial_{t}W)(\partial_{t}{\cal Y})\sim\mathcal{O}(\xi^{3})$ .

Regarding the contribution from the shear stress, we simply follow the Schumaker-Thorne formalism [6], which employs the Hookean approximation so that the shear contribution to the stress-energy tensor from the shear tensor, $S_{\mu\nu}$ , is proportional only to the shear modulus, i.e.,

T_{\mu\nu}^{\rm(S)}=-2\mu S_{\mu\nu}.

(10)

The shear tensor is related to the “rate of shear” tensor, $\sigma_{\mu\nu}$ , [56]

{\cal L}_{u}S_{\mu\nu}=\frac{2}{3}S_{\mu\nu}\nabla_{\alpha}u^{\alpha}+\sigma_{\mu\nu},

(11)

where ${\cal L}_{u}$ denotes the Lie derivative along the direction of four-velocity, $u^{\mu}$ , and one has

	$\displaystyle{\cal L}_{u}S_{\mu\nu}$	$\displaystyle=$	$\displaystyle u^{\alpha}\nabla_{\alpha}S_{\mu\nu}+S_{\alpha\nu}\nabla_{\mu}u^{\alpha}+S_{\mu\alpha}\nabla_{\nu}u^{\alpha},$		(12)
	$\displaystyle\sigma_{\mu\nu}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(P^{\alpha}_{\ \nu}\nabla_{\alpha}u_{\mu}+P^{\alpha}_{\ \mu}\nabla_{\alpha}u_{\nu}\right)-\frac{1}{3}P_{\mu\nu}\nabla_{\alpha}u^{\alpha}\,.$		(13)

Here, $P_{\mu\nu}$ is the projection tensor given by

P_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}\,.

(14)

From Eq. (11), one can get

\delta({\cal L}_{u}S_{\mu\nu})=\delta\sigma_{\mu\nu}+\frac{2}{3}\delta S_{\mu\nu}\delta u^{\alpha}_{\ ;\alpha},

(15)

because $u^{\alpha}_{\ ;\alpha}=0$ and $S_{\mu\nu}=0$ for the (background) equilibrium model. Running through the components $(\mu,\nu)=(\mu,\phi)$ in Eq. (15), it is possible to find

$\displaystyle\delta\sigma_{t\phi}+\frac{2}{3}\delta S_{t\phi}\delta u^{\alpha}_{\ ;\alpha}=$	$\displaystyle e^{-\Phi}\delta S_{t\phi,t}+\delta u^{r}\delta S_{t\phi,r}+\delta u^{\theta}\delta S_{t\phi,\theta}+\delta S_{r\phi}\delta u^{r}_{\ ,t}+\delta S_{\theta\phi}\delta u^{\theta}_{\ ,t}+\delta S_{\phi\phi}\delta u^{\phi}_{\ ,t},$	(16)
$\displaystyle\delta\sigma_{r\phi}+\frac{2}{3}\delta S_{r\phi}\delta u^{\alpha}_{\ ;\alpha}=$	$\displaystyle e^{-\Phi}(\delta S_{r\phi,t}-\Phi_{,r}\delta S_{t\phi})+\delta u^{r}\delta S_{r\phi,r}+\delta u^{\theta}\delta S_{r\phi,\theta}+\delta S_{r\phi}\delta u^{r}_{\ ,r}+\delta S_{\theta\phi}\delta u^{\theta}_{\ ,r}+\delta S_{\phi\phi}\delta u^{\phi}_{\ ,r},$	(17)
$\displaystyle\delta\sigma_{\theta\phi}+\frac{2}{3}\delta S_{\theta\phi}\delta u^{\alpha}_{\ ;\alpha}=$	$\displaystyle e^{-\Phi}\delta S_{\theta\phi,t}+\delta u^{r}\delta S_{\theta\phi,r}+\delta u^{\theta}\delta S_{\theta\phi,\theta}+\delta S_{r\phi}\delta u^{r}_{\ ,\theta}+\delta S_{\theta\phi}\delta u^{\theta}_{\ ,\theta}+\delta S_{\phi\phi}\delta u^{\phi}_{\ ,\theta},$	(18)
$\displaystyle\delta\sigma_{\phi\phi}+\frac{2}{3}\delta S_{\phi\phi}\delta u^{\alpha}_{\ ;\alpha}=$	$\displaystyle e^{-\Phi}\delta S_{\phi\phi,t}+\delta u^{r}\delta S_{\phi\phi,r}+\delta u^{\theta}\delta S_{\phi\phi,\theta},$	(19)

where $\delta\sigma_{\mu\phi}$ and $\delta u^{\alpha}_{\ ;\alpha}$ up to 2nd order in perturbations are given from the definitions as

$\displaystyle\delta\sigma_{t\phi}=$	$\displaystyle\frac{e^{\Phi}}{2}\big{[}-\Phi_{,r}\delta u^{r}\delta u_{\phi}-\delta u^{r}\delta u_{\phi,r}-\delta u^{\theta}\delta u_{\phi,\theta}\big{]}+\frac{1}{3}e^{\Phi}\delta u_{\phi}\delta u^{\alpha}_{\ ;\alpha},$	(20)
$\displaystyle\delta\sigma_{r\phi}=$	$\displaystyle\frac{1}{2}\left[\Phi_{,r}\delta u_{\phi}-\frac{2}{r}\delta u_{\phi}+\delta u_{\phi,r}+e^{-\Phi}\delta u_{\phi}\delta u_{r,t}+e^{-\Phi}\delta u_{r}\delta u_{\phi,t}\right],$	(21)
$\displaystyle\delta\sigma_{\theta\phi}=$	$\displaystyle\frac{1}{2}\left[\delta u_{\phi,\theta}-2\frac{\cos\theta}{\sin\theta}\delta u_{\phi}+e^{-\Phi}\delta u_{\phi}\delta u_{\theta,t}+e^{-\Phi}\delta u_{\theta}\delta u_{\phi,t}\right],$	(22)
$\displaystyle\delta\sigma_{\phi\phi}=$	$\displaystyle re^{-2\Lambda}\sin^{2}\theta\delta u_{r}+\sin\theta\cos\theta\delta u_{\theta}-\frac{1}{3}r^{2}\sin^{2}\theta\delta u^{\alpha}_{\ ;\alpha}+e^{-\Phi}\delta u_{\phi}\delta u_{\phi,t},$	(23)
$\displaystyle\delta u^{\alpha}_{\ ;\alpha}=$	$\displaystyle re^{-\Phi}\partial_{tr}W+re^{-\Phi}\left(\Lambda_{,r}+\frac{3}{r}\right)\partial_{t}W+e^{-\Phi}\partial_{t\theta}V+\frac{\cos\theta}{\sin\theta}e^{-\Phi}\partial_{t}V.$	(24)

III.1 Order matching

At this stage, it is instructive to expand the relevant quantities up to successive orders. We recall that we consider only axial perturbations up to second-order, with the polar parity ones truncated at linear order. That is, $W^{(2)}=V^{(2)}=0$ but ${\cal Y}^{(2)}\neq 0$ . Some algebra now reveals that

$\displaystyle\delta\sigma_{t\phi}^{(1)}=$	$\displaystyle 0\,,$	(25)
$\displaystyle\delta\sigma_{r\phi}^{(1)}=$	$\displaystyle\frac{r^{2}}{2}e^{-\Phi}\sin\theta(\partial_{tr}{\cal Y}^{(1)})\,,$	(26)
$\displaystyle\delta\sigma_{\theta\phi}^{(1)}=$	$\displaystyle\frac{r^{2}}{2}e^{-\Phi}\left[\sin\theta(\partial_{t\theta}{\cal Y}^{(1)})-\cos\theta(\partial_{t}{\cal Y}^{(1)})\right]\,,$	(27)
$\displaystyle\delta\sigma_{\phi\phi}^{(1)}=$	$\displaystyle-\frac{r^{2}}{3}e^{-\Phi}\sin^{2}\theta\left[r\partial_{tr}W^{(1)}+r\Lambda_{,r}\partial_{t}W^{(1)}+\partial_{t\theta}V^{(1)}-\frac{2\cos\theta}{\sin\theta}\partial_{t}V^{(1)}\right]\,,$	(28)
$\displaystyle\delta u^{\alpha(1)}_{\ ;\alpha}=$	$\displaystyle re^{-\Phi}\partial_{tr}W^{(1)}+re^{-\Phi}\left(\Lambda_{,r}+\frac{3}{r}\right)\partial_{t}W^{(1)}+e^{-\Phi}\partial_{t\theta}V^{(1)}+\frac{\cos\theta}{\sin\theta}e^{-\Phi}\partial_{t}V^{(1)}\,,$	(29)

while

$\displaystyle\delta\sigma_{t\phi}^{(2)}=$	$\displaystyle{\mathfrak{s}}_{0}({\cal P},{\cal A})\,,$	(30)
$\displaystyle\delta\sigma_{r\phi}^{(2)}=$	$\displaystyle\frac{r^{2}}{2}e^{-\Phi}\sin\theta(\partial_{tr}{\cal Y}^{(2)})+{\mathfrak{s}}_{1}({\cal P},{\cal A})\,,$	(31)
$\displaystyle\delta\sigma_{\theta\phi}^{(2)}=$	$\displaystyle\frac{r^{2}}{2}e^{-\Phi}\left[\sin\theta(\partial_{t\theta}{\cal Y}^{(2)})-\cos\theta(\partial_{t}{\cal Y}^{(2)})\right]+{\mathfrak{s}}_{2}({\cal P},{\cal A})\,,$	(32)
$\displaystyle\delta\sigma_{\phi\phi}^{(2)}=$	$\displaystyle re^{-2\Lambda}\sin^{2}\theta\delta u_{r}^{(2)}+\sin\theta\cos\theta\delta u_{\theta}^{(2)}-\frac{1}{3}r^{2}\sin^{2}\theta\delta u^{\alpha(2)}_{\ ;\alpha}+e^{-\Phi}\delta u_{\phi}^{(1)}\delta u_{\phi,t}^{(1)}\,,$	(33)

where ${\mathfrak{s}}_{0}({\cal P},{\cal A})$ , ${\mathfrak{s}}_{1}({\cal P},{\cal A})$ , and ${\mathfrak{s}}_{2}({\cal P},{\cal A})$ are the terms composed of the combination of the linear polar ( $\mathcal{P}=V^{(1)},W^{(1)}$ ) and linear axial perturbations ( $\mathcal{A}=\mathcal{Y}^{(1)}$ ), whose forms are concretely given in Appendix A.

At linear-order, Eqs. (16) – (19) become

$\displaystyle\delta S^{(1)}_{t\phi}$	$\displaystyle=$	$\displaystyle 0\,,$	(34)
$\displaystyle\delta S^{(1)}_{r\phi}$	$\displaystyle=$	$\displaystyle\frac{r^{2}}{2}\sin\theta(\partial_{r}{\cal Y}^{(1)})\,,$	(35)
$\displaystyle\delta S^{(1)}_{\theta\phi}$	$\displaystyle=$	$\displaystyle\frac{r^{2}}{2}\left[\sin\theta(\partial_{\theta}{\cal Y}^{(1)})-\cos\theta\,{\cal Y}^{(1)}\right]\,,$	(36)
$\displaystyle\delta S^{(1)}_{\phi\phi}$	$\displaystyle=$	$\displaystyle-\frac{r^{2}}{3}\sin^{2}\theta\left[r\partial_{r}W^{(1)}+r\Lambda_{,r}W^{(1)}+\partial_{\theta}V^{(1)}-\frac{2\cos\theta}{\sin\theta}V^{(1)}\right].$	(37)

We assume that $\delta S^{(1)}_{\mu\nu}=0$ at $t=0$ , i.e., stress builds up from zero from a relaxed background state. We note that $\delta S_{r\phi}^{(1)}$ and $\delta S_{\theta\phi}^{(1)}$ are only composed of linear axial perturbations, while $\delta S_{\phi\phi}^{(1)}$ is only composed of linear polar perturbations. Further note that the other components $\delta S^{(1)}_{ij}$ are non-zero, but do not contribute to the torsional-mode dynamics in a static, unmagnetized star; they can be found elsewhere in the literature (e.g. [26]).

Moving on to second-order, Eqs. (16) – (18) become

$\displaystyle\delta S_{t\phi,t}^{(2)}=$	$\displaystyle e^{\Phi}\left[\delta\sigma_{t\phi}^{(2)}-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,t}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,t}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,t}\right],$	(38)
$\displaystyle\delta S_{r\phi,t}^{(2)}=$	$\displaystyle\Phi_{,r}\delta S_{t\phi}^{(2)}+e^{\Phi}\left[\delta\sigma_{r\phi}^{(2)}+\frac{2}{3}\delta S_{r\phi}^{(1)}\delta u^{\alpha(1)}_{\ ;\alpha}-\delta u^{r(1)}\delta S_{r\phi,r}^{(1)}-\delta u^{\theta(1)}\delta S_{r\phi,\theta}^{(1)}-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,r}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,r}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,r}\right],$	(39)
$\displaystyle\delta S_{\theta\phi,t}^{(2)}=$	$\displaystyle e^{\Phi}\left[\delta\sigma_{\theta\phi}^{(2)}+\frac{2}{3}\delta S_{\theta\phi}^{(1)}\delta u^{\alpha(1)}_{\ ;\alpha}-\delta u^{r(1)}\delta S_{\theta\phi,r}^{(1)}-\delta u^{\theta(1)}\delta S_{\theta\phi,\theta}^{(1)}-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,\theta}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,\theta}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,\theta}\right].$	(40)

We note that the second-order perturbation, $\delta S_{\phi\phi}^{(2)}$ , does not contribute to the second-order equations for torsional oscillations. From Eqs. (38) – (40), one can derive

$\displaystyle\delta S_{t\phi,t}^{(2)}\equiv$	$\displaystyle{\mathfrak{S}}_{0}({\cal P},{\cal A}),$	(41)
$\displaystyle\delta S_{r\phi,t}^{(2)}=$	$\displaystyle\frac{r^{2}}{2}\sin\theta(\partial_{tr}{\cal Y}^{(2)})+{\mathfrak{S}}_{1}({\cal P},{\cal A})+\Phi_{,r}\delta S_{t\phi}^{(2)},$	(42)
$\displaystyle\delta S_{\theta\phi,t}^{(2)}=$	$\displaystyle\frac{r^{2}}{2}\left[\sin\theta(\partial_{t\theta}{\cal Y}^{(2)})-\cos\theta(\partial_{t}{\cal Y}^{(2)})\right]+{\mathfrak{S}}_{2}({\cal P},{\cal A}),$	(43)

where again the Gothic symbols ${\mathfrak{S}}_{0}({\cal P},{\cal A})$ , ${\mathfrak{S}}_{1}({\cal P},{\cal A})$ , and ${\mathfrak{S}}_{2}({\cal P},{\cal A})$ represent terms composed of combinations of the linear polar and linear axial perturbations; they are shown explicitly in Appendix A. These equations can be immediately integrated, again imposing the initial condition that we start with an unstressed state at all orders. By defining ${{\cal S}}_{i}({\cal P},{\cal A})$ for $i=0,1,2$ as

\displaystyle{\cal S}_{i}({\cal P},{\cal A})\equiv\int_{0}^{t}{\mathfrak{S}}_{i}({\cal P},{\cal A})dt,

(44)

one can get

$\displaystyle\delta S^{(2)}_{t\phi}=$	$\displaystyle{{\cal S}}_{0}({\cal P},{\cal A})\,,$	(45)
$\displaystyle\delta S^{(2)}_{r\phi}=$	$\displaystyle\frac{r^{2}}{2}\sin\theta(\partial_{r}{\cal Y}^{(2)})+{\cal S}_{1}({\cal P},{\cal A})+\Phi_{,r}\int_{0}^{t}\delta S_{t\phi}^{(2)}dt\,,$	(46)
$\displaystyle\delta S^{(2)}_{\theta\phi}=$	$\displaystyle\frac{r^{2}}{2}\left[\sin\theta(\partial_{\theta}{\cal Y}^{(2)})-\cos\theta\,{\cal Y}^{(2)}\right]+{\cal S}_{2}({\cal P},{\cal A})\,.$	(47)

Now, simply⁴⁴4In considering second-order perturbations of the shear portion of the energy-momentum tensor we could, in principle, have terms that are proportional to $\delta\mu$ . Such terms could be self-consistently associated with the density perturbation, $\delta\varepsilon$ , implicitly through Eq. (3). This is not trivial however since the shear modulus depends only on the ion density, rather than the overall density, $\delta\varepsilon$ . In practice, if one considers a contribution from $\delta\mu$ , one has to add the terms $\delta\mu^{(1)}\delta S_{r\phi}^{(1)}$ to Eq. (52) and $\delta\mu^{(1)}\delta S_{\theta\phi}^{(1)}$ to Eq. (53). These additional terms change the source term ${\cal S}_{S}({\cal P},{\cal A})$ in Eq. (58) a little, but as mentioned below, the main results shown in this study are unchanged. considering $\delta T^{\rm(S)}_{\mu\nu}=-2\mu\delta S_{\mu\nu}$ from Eq. (10), i.e., assuming that $\delta\mu=0$ , we find its linear components

$\displaystyle\delta T^{{\rm(S)}(1)}_{t\phi}=$	$\displaystyle 0\,,$	(48)
$\displaystyle\delta T^{{\rm(S)}(1)}_{r\phi}=$	$\displaystyle-\mu r^{2}\sin\theta(\partial_{r}{\cal Y}^{(1)})\,,$	(49)
$\displaystyle\delta T^{{\rm(S)}(1)}_{\theta\phi}=$	$\displaystyle-\mu r^{2}\left[\sin\theta(\partial_{\theta}{\cal Y}^{(1)})-\cos\theta\,{\cal Y}^{(1)}\right]\,,$	(50)

as well as its second-order components

$\displaystyle\delta T^{{\rm(S)}(2)}_{t\phi}=$	$\displaystyle-2\mu{\cal S}_{0}({\cal P},{\cal A})\,,$	(51)
$\displaystyle\delta T^{{\rm(S)}(2)}_{r\phi}=$	$\displaystyle-\mu r^{2}\sin\theta(\partial_{r}{\cal Y}^{(2)})-2\mu{\cal S}_{1}({\cal P},{\cal A})-2\mu\Phi_{,r}\int_{0}^{t}\delta S_{t\phi}^{(2)}dt,$	(52)
$\displaystyle\delta T^{{\rm(S)}(2)}_{\theta\phi}=$	$\displaystyle-\mu r^{2}\left[\sin\theta(\partial_{\theta}{\cal Y}^{(2)})-\cos\theta\,{\cal Y}^{(2)}\right]-2\mu{\cal S}_{2}({\cal P},{\cal A})\,.$	(53)

The divergence of $\delta T^{{\rm(S)}(1)\phi\nu}$ is given as

	$\displaystyle\delta T^{{\rm(S)}(1)\phi\nu}_{\qquad\quad;\nu}=$	$\displaystyle-\left[\frac{\mu e^{-2\Lambda}}{\sin\theta}(\partial_{r}{\cal Y}^{(1)})\right]_{,r}+\frac{\mu}{r^{2}}\left[\frac{\cos\theta}{\sin^{2}\theta}{\cal Y}^{(1)}-\frac{1}{\sin\theta}(\partial_{\theta}{\cal Y}^{(1)})\right]_{,\theta}$
		$\displaystyle-\frac{\mu e^{-2\Lambda}}{\sin\theta}\left(\Phi_{,r}+\Lambda_{,r}+\frac{4}{r}\right)(\partial_{r}{\cal Y}^{(1)})+\frac{3\mu}{r^{2}}\left[\frac{\cos^{2}\theta}{\sin^{3}\theta}{\cal Y}^{(1)}-\frac{\cos\theta}{\sin^{2}\theta}(\partial_{\theta}{\cal Y}^{(1)})\right]\,.$		(54)

As is well-known, the relevant equations of motion are separable at first-order. Abusing notation slightly and decomposing ${\cal Y}^{(1)}(t,r,\theta)=\sum_{\ell}{\cal Y}_{\ell}^{(1)}(t,r)\partial_{\theta}P_{\ell}(\cos\theta)$ for Legendre polynomials $P_{\ell}(\cos\theta)$ , one can express the divergence $\delta T^{{\rm(S)}(1)\phi\nu}_{\ \ \ \ \ \ \ \ ;\nu}$ , for each $\ell$ , as

\displaystyle\delta T^{{\rm(S)}(1)\phi\nu}_{\ \ \ \ \ \ \ \ ;\nu}=

\displaystyle-\mu e^{-2\Lambda}\left[\partial_{rr}{\cal Y}_{\ell}^{(1)}+\left(\frac{\mu_{,r}}{\mu}+\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)(\partial_{r}{\cal Y}_{\ell}^{(1)})-\frac{(\ell+2)(\ell-1)}{r^{2}}e^{2\Lambda}{\cal Y}_{\ell}^{(1)}\right]\frac{1}{\sin\theta}\partial_{\theta}P_{\ell},

(55)

where we have used the well known relation $\left[\partial^{2}_{\theta}+(\cos\theta/\sin\theta)\partial_{\theta}+\ell(\ell+1)\right]P_{\ell}=0$ . Thus, the perturbation equation at linear level reduces to

$\displaystyle 0=$	$\displaystyle\delta T^{{\rm(F)}(1)\phi\nu}_{\ \ \ \ \ \ \ \ ;\nu}+\delta T^{{\rm(S)}(1)\phi\nu}_{\ \ \ \ \ \ \ \ ;\nu}$
$\displaystyle=$	$\displaystyle(\varepsilon+p)e^{-2\Phi}(\partial_{tt}{\cal Y}_{\ell}^{(1)})\frac{1}{\sin\theta}\partial_{\theta}P_{\ell}$
	$\displaystyle-\mu e^{-2\Lambda}\left[\partial_{rr}{\cal Y}_{\ell}^{(1)}+\left(\frac{\mu_{,r}}{\mu}+\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)(\partial_{r}{\cal Y}_{\ell}^{(1)})-\frac{(\ell+2)(\ell-1)}{r^{2}}e^{2\Lambda}{\cal Y}_{\ell}^{(1)}\right]\frac{1}{\sin\theta}\partial_{\theta}P_{\ell},$	(56)

and hence

\displaystyle 0=-\frac{\varepsilon+p}{\mu}e^{-2\Phi+2\Lambda}(\partial_{tt}{\cal Y}_{\ell}^{(1)})+\partial_{rr}{\cal Y}_{\ell}^{(1)}+\left(\frac{\mu_{,r}}{\mu}+\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)(\partial_{r}{\cal Y}_{\ell}^{(1)})-\frac{(\ell+2)(\ell-1)}{r^{2}}e^{2\Lambda}{\cal Y}_{\ell}^{(1)},

(57)

as in Ref. [6].

On the other hand, the divergence of $\delta T^{{\rm(S)}\phi\nu}$ at second-order reads

	$\displaystyle\delta T^{{\rm(S)}(2)\phi\nu}_{\ \ \ \ \ \ \ \ ;\nu}=$	$\displaystyle\frac{\mu}{r^{2}\sin\theta}\bigg{[}-r^{2}e^{-2\Lambda}(\partial_{rr}{\cal Y}^{(2)})-r^{2}e^{-2\Lambda}\left(\frac{\mu_{,r}}{\mu}+\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)(\partial_{r}{\cal Y}^{(2)})$
	$\displaystyle-$	$\displaystyle\partial_{\theta\theta}{\cal Y}^{(2)}-\cot\theta\partial_{\theta}{\cal Y}^{(2)}+\left(\cot^{2}\theta-1\right){\cal Y}^{(2)}\bigg{]}+{\cal S}_{S}({\cal P},{\cal A})\,,$		(58)

where ${\cal S}_{S}({\cal P},{\cal A})$ is a coupling term composed of combinations of linear polar and axial perturbations, whose form is explicitly shown in Appendix A. Thus, the second-order perturbation equation for axial-parity oscillations can be derived from

\displaystyle\delta T^{(2)\phi\nu}_{\qquad;\nu}=\delta T^{{\rm(F)}(2)\phi\nu}_{\qquad\,;\nu}+\delta T^{{\rm(S)}(2)\phi\nu}_{\qquad\,;\nu}=0,

(59)

which leads to

	$\displaystyle-(\varepsilon+p)e^{-2\Phi+2\Lambda}(\partial_{tt}{\cal Y}^{(2)})+$	$\displaystyle\mu(\partial_{rr}{\cal Y}^{(2)})+\left[\mu_{,r}+\mu\left(\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)\right](\partial_{r}{\cal Y}^{(2)})$
	$\displaystyle+$	$\displaystyle\frac{\mu}{r^{2}}e^{2\Lambda}\left[\partial_{\theta\theta}{\cal Y}^{(2)}+\cot\theta\partial_{\theta}{\cal Y}^{(2)}+\left(1-\cot^{2}\theta\right){\cal Y}^{(2)}\right]=e^{2\Lambda}\sin\theta\left[{\cal S}_{F}({\cal P},{\cal A})+{\cal S}_{S}({\cal P},{\cal A})\right],$		(60)

where ${\cal S}_{F}({\cal P},{\cal A})$ is another part of the total coupling term composed of combinations of linear polar and axial perturbations, whose form is explicitly shown in Appendix A. Now, we introduce a new variable, $\tilde{\cal Y}(t,r,\theta)\equiv{\cal Y}(t,r,\theta)/\sin\theta$ , which leads to $\tilde{\cal Y}^{(2)}={\cal Y}^{(2)}/\sin\theta$ . Then, Eq. (60) becomes

	$\displaystyle-(\varepsilon+p)e^{-2\Phi+2\Lambda}(\partial_{tt}\tilde{\cal Y}^{(2)})+$	$\displaystyle\mu(\partial_{rr}\tilde{\cal Y}^{(2)})+\left[\mu_{,r}+\mu\left(\Phi_{,r}-\Lambda_{,r}+\frac{4}{r}\right)\right](\partial_{r}\tilde{\cal Y}^{(2)})$
	$\displaystyle+$	$\displaystyle\frac{\mu}{r^{2}}e^{2\Lambda}\left[(\partial_{\theta\theta}\tilde{\cal Y}^{(2)})+3\partial_{\theta}\tilde{\cal Y}^{(2)}{\cot\theta}\right]=e^{2\Lambda}\left[{\cal S}_{F}({\cal P},{\cal A})+{\cal S}_{S}({\cal P},{\cal A})\right].$		(61)

We note that the operator defining Eq. (61) [i.e. the homogenous equation neglecting source terms] is identical to the linear equation for torsional oscillations. Once one selects specific linear oscillation modes for the axial and polar parity to be excited, one can examine the second-order spectrum. Furthermore, due to the coupling terms, Eq. (61) is generally not separable.

We close this section by stating that, in what remains of this paper, we numerically evolve Eq. (60) in isolation. That is, we do not solve a closed system where mode amplitudes adjust iteratively as seed energy is sapped from the first-order spectrum to grow modes at second-order as has been done, for example, in the $f$ - [59, 60] and $r$ -mode [61, 62] literatures. In reality, Eq. (57) must be solved simultaneously to study the energy-transfer problem and onset of resonant or nonresonant parametric instabilities (i.e., parent-daughter-daughter couplings in addition to the direct parent-parent-daughter couplings). However, such complexity lies beyond the scope of this work, which is mostly concerned with deriving the relevant equations and illustrating how axial-polar couplings at first order can enrich the nonlinear spectrum.

III.2 Boundary and initial conditions

The linear perturbation equation for torsional oscillations, Eq. (57), can be solved as an eigenvalue problem by taking a harmonic time-dependence $e^{i\omega t}$ for angular frequency $\omega=2\pi f$ . We impose a zero traction condition at the base of the crust and the zero-torque condition at the interface between the crust and envelope (see e.g., Refs. [6, 24] for details). On the other hand, the equations and boundary conditions for linear polar-parity oscillations are explicitly shown in Ref. [57, 58]. With crustal elasticity, the interface ( $i_{i}$ -) and shear ( $s_{i}$ -) modes can be excited in addition to the usual acoustic oscillations, i.e., the fundamental ( $f$ -) and pressure ( $p_{i}$ -) modes. Since the number of excitable $i$ -modes is generally the same as the number of interfaces, where the non-zero shear modulus discontinuously drops to zero, we can expect that two $i$ -modes exist in our case. In addition, since we simply consider linear polar perturbations without stratification for zero-temperature nuclear matter in this study, we cannot discuss couplings induced by the gravity ( $g$ -) modes. However, high-overtones of $g$ -modes are of low frequency, which may play an important role in explaining the observed QPOs in magnetar flares. Their inclusion is relatively straightforward, and will be investigated elsewhere.

Since it is not a priori clear which multipoles $\ell$ will be excited by non-linear effects for some given initial data, we examine oscillations inside the crustal region over the whole angular domain $0\leq\theta\leq\pi$ in this study. To evolve Eq. (61) in two-dimensional space ( $r$ , $\theta$ ), one has to impose some boundary conditions. Here, as in the case of the linear perturbations with axial parity, we impose a zero traction condition at the base of the crust ( $r=R_{c}$ ) and at the interface between the crust and envelope ( $r=R_{e}$ ). In symbols, we enforce

\partial_{r}\tilde{\cal Y}^{(2)}|_{r=R_{c},R_{e}}=0,

(62)

because of Eq. (52). On the other hand, since we simply consider axisymmetric perturbations, the boundary conditions imposed at $\theta=0$ and $\pi$ should be $\partial_{\theta}\tilde{\cal Y}^{(2)}|_{\theta=0,\pi}=0$ . The time evolution is calculated with the iterated Crank-Nicholson method [63]. Hereafter, $N_{r}$ and $N_{\theta}$ denote the grid number in the radial and $\theta$ direction, respectively.

In addition, we simply assume that ${\cal Y}^{(2)}=0$ at $t=0$ as an initial condition. That is, each ${\cal Y}^{(2)}$ solution shown below is completely induced by the source terms composed of the product of the linear axial and polar perturbations, with the coupling timescale being essentially instantaneous.

IV Results: Nonlinear Coupling

The numerical results strongly depend on the spatial resolution. Numerical convergence tests, carried out by evolving Eq. (61) assuming that ${\cal S}_{F}({\cal P},{\cal A})={\cal S}_{S}({\cal P},{\cal A})=0$ (i.e., the linear problem), are presented in Appendix B. These demonstrate how the results vary with radial, $N_{r}$ , and angular, $N_{\theta}$ , resolutions. Based on that analysis, we find that one adequately resolves the modes up to the 1st overtones with angular quantum-numbers reaching $\ell\simeq 10$ for $N_{r}\simeq 500$ and $N_{\theta}\simeq 100$ . So, hereafter we adopt $N_{r}=500$ and $N_{\theta}=120$ for evolutions with nonlinear coupling. We note that higher-order ( $n>1$ ) overtones may not be resolved with our resolution even if they are excited, but these lie at far-too-high frequencies to be observable. If considering problems with energy-transfer, these could still be relevant however (as for, e.g., $p$ - $g$ couplings [43]).

We now turn to the main results of this work and examine the torsional oscillations induced by non-linear couplings between linear perturbations of both axial and polar parity by evolving Eq. (61) with the source terms, ${\cal S}_{F}({\cal P},{\cal A})$ and ${\cal S}_{S}({\cal P},{\cal A})$ . As mentioned above, order- $\ell$ linear perturbations are independent of all other order- $\ell^{\prime}\neq\ell$ linear perturbations. The polar perturbations are also independent of the axial ones. So, one can examine the pattern of mode excitation induced by the linear coupling via Eq. (61) by selecting a specific combination of axial and polar oscillation modes, e.g., ${}_{2}t_{0}$ (the $\ell=2$ fundamental torsional mode) and ${}_{2}s_{1}$ (the $\ell=2$ first shear mode). In this study, as a first step, we consider the coupling between the torsional oscillations, ${}_{\ell}t_{0}$ (or ${}_{\ell}t_{1}$ ), and the $\ell=2$ polar oscillations (though we could consider modes of multipolarity even higher than $\ell=11$ for $N_{\theta}=120$ , as detailed in Appendix B). A study involving other couplings may be shown elsewhere.

Several eigenfrequencies determined via the eigenvalue problem with a linear perturbation analysis, using the neutron star model constructed in Sec. II, are listed in Table 1. We note that in this study we identify the $i_{1}$ - and $i_{2}$ -modes in such a way that the $i_{2}$ -mode is the mode excited at the crust-envelope interface, while the $i_{1}$ -mode (tangential) eigenfunction has jumps at both the crust-core and crust-envelope interfaces. With the source terms, ${\cal S}_{F}({\cal P},{\cal A})$ and ${\cal S}_{S}({\cal P},{\cal A})$ , determined explicitly using the quoted eigenfrequencies and associated eigenfunctions, we evolve Eq. (61) for 10 seconds. For this, we expect a bandwidth of $\sim 0.1~{}$ Hz. Even if we select a specific combination of linear axial and polar oscillation modes, there is still freedom in setting their amplitudes. In this study, we simply set the maximum amplitude of the selected linear axial ( ${\cal Y}^{(1)}$ ) and polar perturbations (either $W^{(1)}$ or $V^{(1)}$ ) in the elastic region to be unity. On the other hand, since the source terms, ${\cal S}_{F}({\cal P},{\cal A})$ and ${\cal S}_{S}({\cal P},{\cal A})$ , are always composed of the product of the amplitudes of linear axial and polar perturbations, it is not the individual amplitudes that are important but rather their product in the evolution of ${\cal Y}^{(2)}$ . One can expect the results shown in this study would be qualitatively unchanged, even if one changes such a product into some value less than 1.

Table 1: Eigenfrequencies determined via the (first-order) eigenvalue problem for the neutron star model considered in this study, i.e.,

M=1.41M_{\odot}

and

R=11.68

km constructed with the SLy4 EOS. The labelling of these modes is defined in Sec. III.

Mode	Frequency (Hz)
${}_{2}t_{0}$	23.9
${}_{3}t_{0}$	37.8
${}_{4}t_{0}$	50.7
${}_{5}t_{0}$	63.2
${}_{6}t_{0}$	75.6
${}_{2}t_{1}$	898.5
${}_{3}t_{1}$	899.0
${}_{4}t_{1}$	899.7
${}_{5}t_{1}$	900.6
${}_{6}t_{1}$	901.7
${}_{2}i_{2}$	40.9
${}_{2}i_{1}$	45.2
${}_{2}s_{1}$	895.1
${}_{2}s_{2}$	1463.8
${}_{2}s_{3}$	1806.5
${}_{2}s_{4}$	2276.9
${}_{2}s_{5}$	2829.0
${}_{2}f$	2407.1

IV.1 Seeding via coupled lowest- $n$ modes

First, we show the results with coupling between ${}_{2}t_{0}$ and ${}_{2}s_{1}$ . The FFT computed via the resulting waveform for $\mathcal{Y}^{(2)}$ is shown in Fig. 1. We find several modes are excited, which we identify as ${}_{2}t_{0}$ , ${}_{4}t_{0}$ , ${}_{2}t_{1}$ , and ${}_{4}t_{1}$ . In addition to these frequencies, which appear in the linear perturbation analysis, we also find additional modes, whose frequencies correspond to the sum of that associated with ${}_{2}t_{0}$ and ${}_{2}s_{1}$ . We note that the oscillations with ${}_{2}t_{1}$ , ${}_{4}t_{1}$ , and ${}_{2}t_{0}+{}_{2}s_{1}$ (= 919.0 Hz) seem to be excited more strongly than those with ${}_{2}t_{0}$ and ${}_{4}t_{0}$ at the 2nd-order perturbation level, even though we consider the ${}_{2}t_{0}$ mode to build the source terms. Note we are showing only the second-order spectrum, not FFTs of the total eigenfunction $\mathcal{Y}$ . More precisely, although we sometimes find $\mathcal{Y}^{(2)}_{\rm amp,daughter}>\mathcal{Y}^{(2)}_{\rm amp,parents}$ , we still have $\mathcal{Y}^{(1)}_{\rm amp,daughter}+\mathcal{Y}^{(2)}_{\rm amp,daughter}=\mathcal{Y}^{(2)}_{\rm amp,daughter}\ll\mathcal{Y}^{(1)}_{\rm amp,parents}+\mathcal{Y}^{(2)}_{\rm amp,parents}$ , where $\mathcal{Y}^{(i)}_{\rm amp,daughter}$ and $\mathcal{Y}^{(i)}_{\rm amp,parents}$ denote the amplitude of the $i$ -th order perturbations for the daughter and parent modes, respectively.

Refer to caption — Figure 1: The FFT of $\mathcal{Y}^{(2)}$ evolved for 10 seconds. Here we imposed a coupling (i.e. initial data) between ${}_{2}t_{0}$ (23.9 Hz) and ${}_{2}s_{1}$ (895.1 Hz). For reference, the ${}_{2}t_{0}$ and ${}_{2}s_{1}$ frequencies are shown with dotted lines and we have focused on the regions $0\leq f\leq 200$ Hz and $880\leq f\leq 960$ Hz as no modes are excited in the interval $200\leq f\leq 880$ Hz.

The reason why ${}_{3}t_{0}$ is not excited can be understood by carefully examining the perturbation equations (61). In general, the odd- and even-order Legendre polynomials for $\ell\geq 2$ can be written as

	$\displaystyle P_{2m}=\sum_{j=0}^{m}a_{j}\cos^{2j}\theta,$		(63)
	$\displaystyle P_{2m+1}=\sum_{j=0}^{m}b_{j}\cos^{2j+1}\theta,$		(64)

where $a_{j}$ and $b_{j}$ are appropriate constants. Meanwhile, the angular profile of polar and axial perturbations are essentially expressed as $P_{\ell_{p}}$ and $\partial_{\theta}P_{\ell_{a}}/\sin\theta$ , respectively, where the indices $\ell$ for the polar and axial perturbations are explicitly shown as $\ell_{p}$ and $\ell_{a}$ . So, the coupling between the linear polar and linear axial perturbations in the source terms, ${\cal S}_{i=F,S}$ in Eq. (61) for a single mode behave like

$\displaystyle{\cal S}_{i}({\cal P},{\cal A})$	$\displaystyle\simeq P_{\ell_{p}}\frac{\partial_{\theta}P_{\ell_{a}}}{\sin\theta}$	(65)
	$\displaystyle\simeq\sum_{j=0}^{m_{p}+m_{a}}\alpha_{j}\left(\frac{\partial_{\theta}P_{2j}}{\sin\theta}\right)\ \ \ {\rm for}\ (\ell_{p},\ell_{a})=(2m_{p},2m_{a}),$	(66)
	$\displaystyle\simeq\sum_{j=0}^{m_{p}+m_{a}}\alpha_{j}\left(\frac{\partial_{\theta}P_{2j+1}}{\sin\theta}\right)\ \ \ {\rm for}\ (\ell_{p},\ell_{a})=(2m_{p},2m_{a}+1),$	(67)
	$\displaystyle\simeq\sum_{j=0}^{m_{p}+m_{a}}\alpha_{j}\left(\frac{\partial_{\theta}P_{2j+1}}{\sin\theta}\right)\ \ \ {\rm for}\ (\ell_{p},\ell_{a})=(2m_{p}+1,2m_{a}),$	(68)
	$\displaystyle\simeq\sum_{j=0}^{m_{p}+m_{a}+1}\alpha_{j}\left(\frac{\partial_{\theta}P_{2j}}{\sin\theta}\right)\ \ \ {\rm for}\ (\ell_{p},\ell_{a})=(2m_{p}+1,2m_{a}+1),$	(69)

where $\alpha_{j}$ is again some constant which is a combination of $a_{j}$ and $b_{j}$ . That is, the source terms should induce axial-type oscillations up to order $\ell_{p}+\ell_{a}$ . Moreover, from the above expression, we find that torsional oscillations with only even (odd) values of $\ell$ are selectively excited when one considers polar and axial couplings such that $\ell_{p}$ and $\ell_{a}$ are both even or both odd numbers (even and odd numbers or odd and even numbers). For example, only the $\ell=2$ , 4, and 6 torsional oscillations can be excited at second-order by the coupling between $\ell_{p}=2$ and $\ell_{a}=4$ at first order; similarly, only the $\ell=3$ and 5 second-order torsional oscillations can be excited by the coupling between $\ell_{p}=2$ and $\ell_{a}=3$ .

In addition to these mode excitations, since the source term in Eq. (61) is composed of a combination of linear polar and linear axial perturbations, which respectively depend on time as $\exp(i\omega_{p}t)$ and $\exp(i\omega_{a}t)$ with eigenfrequencies $\omega_{p}$ and $\omega_{a}$ of the considered linear perturbations, one can expect an additional mode with frequency $\omega=\omega_{p}+\omega_{a}$ will be excited. This is simply a property of inhomogeneous wave equations with harmonic source terms. Note that, in a full treatment of the non-linear problem, one may expect the negative branch, $\omega=|\omega_{p}-\omega_{a}|$ , to be excited also. This is because the eigenvalue equations depend on $\omega^{2}$ , and hence both positive and negative roots are valid solutions. If we were to include source terms for total modes including conjugates, we would excite all branches $(\omega=|\omega_{p}\pm\omega_{a}|)$ , though we set the amplitudes of the conjugates to zero for simplicity.

In a similar way, in Fig. 2, we show the FFT obtained for an evolution triggered by the coupling between ${}_{2}t_{0}$ and ${}_{2}i_{1}$ oscillations (top panel), and between the ${}_{2}t_{0}$ and ${}_{2}i_{2}$ oscillations (bottom panel). As expected, we observe the additional mode excitation with a frequency of ${}_{2}t_{0}+{}_{2}i_{1}$ ( $\approx$ 69.0 Hz) for the coupling between ${}_{2}t_{0}$ and ${}_{2}i_{1}$ , and with a frequency of ${}_{2}t_{0}+{}_{2}i_{2}$ ( $\approx$ 64.8 Hz) for the coupling between ${}_{2}t_{0}$ and ${}_{2}i_{2}$ , highlighted by arrows in the right panels. However, in this case for the coupling with $i$ -mode oscillations, we observe not only the $\ell=2$ and $\ell=4$ modes but also many other, even-order (in $\ell$ ) modes are significantly excited. This result cannot be immediately understood from the analytic arguments surrounding Eq. (66). The richness of the excited spectrum may come from the fact that the frequency of the additionally “superposed” excitation mode, ${}_{2}t_{0}+{}_{2}i_{1}$ or ${}_{2}t_{0}+{}_{2}i_{2}$ , is higher than ${}_{4}t_{0}$ , and consequently the ${}_{\ell}t_{0}$ -modes for $\ell\geq 6$ are also excited. In practice, in the coupling of the ${}_{2}t_{0}$ with the ${}_{2}i$ -modes, the ${}_{4}t_{0}$ and ${}_{6}t_{0}$ modes, which are the two sides of the frequencies of ${}_{2}t_{0}+{}_{2}i_{1}$ or ${}_{2}t_{0}+{}_{2}i_{2}$ , are strongly excited, compared to the other modes.

We also note that excitations at higher frequencies (overtones) have more power in the FFT for cases coupled to $s$ -modes, while excitations with lower frequencies (fundamental oscillations) become stronger in cases with a coupling to an $i$ -mode. Due to this feature, as shown in Fig. 3, waveforms in cases with a linear-coupling to an $i$ -mode (left panels) are of lower amplitude (note the scaling on the axes) from that for the coupling with $s$ -mode (right panel). We remind the reader that we have set the (initial) amplitudes as equal for the seeding modes, and these have magnitude 1.

IV.2 Seeding via coupled overtones

Next, in Fig. 4, we show which modes are excited through the coupling of ${}_{2}t_{1}$ (1st overtone of the torsional oscillations, rather than the fundamental mode) with an $i$ -mode [ ${}_{2}i_{2}$ (black) or ${}_{2}i_{1}$ (red)] in the top panel and with an $s$ -mode [ ${}_{2}s_{1}$ (black) or ${}_{2}s_{2}$ (red)] in the bottom panel, focusing on the frequency range lower than 1 kHz.

Focusing first on the upper panel, we see that there are two closeby peaks at $\sim 900$ Hz. One of these corresponds to the ${}_{2}t_{1}$ and ${}_{4}t_{1}$ oscillations ( $\approx 900$ Hz), while the additional frequency denoted with an arrow ( $\approx 950$ Hz), which does not appear in the linear analysis, corresponds to ${}_{2}t_{1}+{}_{2}i_{2}$ ( ${}_{2}t_{1}+{}_{2}i_{1}$ ) for the coupling between ${}_{2}t_{1}$ and ${}_{2}i_{2}$ ( ${}_{2}t_{1}$ and ${}_{2}i_{1}$ ). The former of these persists in the case with a shear mode coupling, while the latter does not. We find that, unlike the case of the coupling of the fundamental torsional oscillations (Fig. 2) with the $i$ -modes, the torsional-mode overtones actually become stronger than the fundamental torsional modes, assuming that the amplitudes of the linear axial and polar perturbations in the source terms are equal to one. This is interesting astrophysically as it suggests that higher-frequency modes might get excited to larger amplitudes when considering higher-frequency seeds. In addition, since the frequency of the additional excitation ( ${}_{2}t_{1}+{}_{2}i_{2}$ or ${}_{2}t_{1}+{}_{2}i_{1}$ ) becomes much higher than those of the fundamental torsional oscillations, the strongest signals in the FFT among the fundamental modes are associated with ${}_{2}t_{0}$ and ${}_{4}t_{0}$ . The feature in the coupling of ${}_{2}t_{1}$ with the $s$ -modes is almost the same as in the case of the coupling of ${}_{2}t_{0}$ , even though the frequency of the additional excitation lies outside of the frequency range considered in the figure.

To consider such phenomena more generally, we show in Fig. 5 the FFTs computed from the waveforms of ${\cal Y}^{(2)}$ when considering seedings from the coupling of ${}_{\ell}t_{0}$ for $\ell=3-6$ with the ${}_{2}i_{2}$ - (top panel), ${}_{2}i_{2}$ - (middle panel), and ${}_{2}s_{i}$ -modes for $i=1,2$ (bottom panel), focusing on the frequency range lower than 200 Hz to best illustrate the spectral structure. One observes a lot of mode excitations as well as the additional mode described earlier (with a frequency of ${}_{\ell}t_{0}+{}_{2}i_{i}$ for $i=1,2$ in the case of a coupling between ${}_{\ell}t_{0}$ and ${}_{2}i$ -modes – couplings with the shear modes are higher frequency). The strongest peaks correspond to this additional mode and those nearest neighbours in terms of frequency. In particular, the ${}_{2}t_{0}$ -mode is not excited for the case of the coupling between ${}_{6}t_{0}$ and ${}_{2}i$ -modes, since the ${}_{2}t_{0}$ -mode frequency is located too far from the additional mode with the frequency of ${}_{6}t_{0}+{}_{2}i_{i}$ for $i=1,2$ . On the other hand, for the case of the coupling between the ${}_{\ell}t_{0}$ and ${}_{2}s$ -modes, one can observe the excitations of the ${}_{\ell^{\prime}}t_{0}$ -modes with up to $\ell^{\prime}=\ell+2$ , as discuss with Eq. (65). Even so, we find that three ${}_{\ell^{\prime}}t_{0}$ -modes with $\ell^{\prime}=\ell-2$ , $\ell$ , and $\ell+2$ become much stronger signals for $\ell\geq 4$ .

On the other hand, Fig. 6 is the same as Fig. 5, but with the coupling of ${}_{\ell}t_{1}$ instead of ${}_{\ell}t_{0}$ for $\ell=3$ , 4, 5 and 6. In this figure, we focus only on the frequency range lower than 200 Hz, but we find that the overtones (and the additional excitation, which becomes in the frequency range of overtone now) are stronger than the fundamental oscillations, as shown in Fig. 4. In the coupling of ${}_{3}t_{1}$ with the ${}_{2}i$ -modes, we find that the modes of ${}_{3}t_{0}$ and ${}_{5}t_{0}$ are strongly excited as expected with Eq. (65), unlike the case for the coupling with the fundamental torsional oscillations. In addition, for $\ell\geq 4$ we find that the modes of ${}_{\ell^{\prime}}t_{0}$ with $\ell^{\prime}=\ell\pm 2$ become strong signals, while the mode with ${}_{\ell}t_{0}$ becomes weak. Meanwhile, in the coupling of ${}_{\ell}t_{1}$ with the ${}_{2}s$ -modes, the feature is more or less similar to the case for the coupling of ${}_{\ell}t_{0}$ as shown in the bottom panel of Fig. 5. We note that the frequency of ${}_{\ell}t_{1}$ weakly depends on the value of $\ell$ , but which modes are excited due to the mode coupling discussed in this study strongly depends on the value of $\ell$ .

V Conclusions

In this paper, we have derived equations describing the evolution of second-order, axisymmetric torsional oscillations inside the elastic crust of a neutron star in general relativity. The relevant master equation, Eq. (60), was solved in the time domain to determine the second-order displacement, ${\cal Y}^{(2)}$ , for a variety of astrophysically-motivated initial conditions. Though we only evolve the second-order functions in the axial sector (i.e. we set $V^{(2)}=W^{(2)}=0$ ), we self-consistently allow ${\cal Y}^{(2)}$ to be sourced by coupled axial-polar oscillations at first order. This represents a useful step forward towards the difficult inverse problem of determining neutron star properties from oscillations observed in the tails of giant flares or other transients. Indeed, most previous studies have focussed on the linear problem, while those that have been carried out at a non-linear level have not included seeding by coupled axial-polar parents. Though we have not considered magnetic fields or crust-core coupling in this paper, already the nonlinear analysis in our relatively simple case reveals some interesting features that are likely to apply to a real, astrophysical star.

One of these features concerns the nature of excited daughters from a given pair of parents. In agreement with the analytic predictions based on Legendre orthogonality (see Eqs. 65–69), a variety of modes can be excited from mixed couplings. Figs. 5 and 6 demonstrate a rich spectrum from just two seed modes of different multipolarity or tonality (torsional plus either interface or shear modes), the amplitudes of which vary by several orders of magnitude and do not necessarily decrease with increasing $\ell$ or $n$ as one might naively expect. This implies that, in a real neutron star system where an elastic overstraining occurs, it is critical to account for nonlinear dynamics.

There are several directions that would be worth pursuing in order to improve on this work. One of these concerns is feedback in the coupling. In this study we have determined the linear oscillation spectrum via the eigenvalue problem, the solutions of which then act as fixed (though still time-dependent) source terms in the second-order equation. However, a more realistic approach would be to loop these into each other to study energy transfer, as in the Newtonian scheme devised by Dziembowski [64] and others (e.g., [59, 60]). By including the energy transfer dynamics one can investigate the saturation amplitude of various excited modes (see also Ref. [65]). These amplitudes can then be compared to the Bayes factors in various statistical analyses for the dynamic spectra in various giant flare tails (e.g. [23]) to approach the inverse problem. Other obvious directions are to include higher-order polar couplings, magnetic fields, and the core. Such studies will be carried out elsewhere.

Acknowledgements.

This work is supported in part by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers JP19KK0354, JP21H01088, JP23K20848, and JP24KF0090, by FY2023 RIKEN Incentive Research Project, and by Pioneering Program of RIKEN for Evolution of Matter in the Universe (r-EMU). This work has been supported by the HORIZON-MSCA-2022 project “EinsteinWaves (101131233)”. AGS is supported by the Prometeo Excellence Programme grant CIPROM/2022/13 from the Generalitat Valenciana.

Appendix A Coupling terms appearing in the main text

In this Appendix, we collectively show the coupling terms appearing in the main text. In the order in which they appear, these read

$\displaystyle{\mathfrak{s}}_{0}({\cal P},{\cal A})=$	$\displaystyle\frac{1}{2}e^{\Phi}\big{[}-\Phi_{,r}\delta u^{r(1)}\delta u_{\phi}^{(1)}-\delta u^{r(1)}\delta u_{\phi,r}^{(1)}-\delta u^{\theta(1)}\delta u_{\phi,\theta}^{(1)}\big{]}+\frac{1}{3}e^{\Phi}\delta u_{\phi}^{(1)}\delta u^{\alpha(1)}_{\ ;\alpha},$	(70)
$\displaystyle{\mathfrak{s}}_{1}({\cal P},{\cal A})=$	$\displaystyle\frac{1}{2}e^{-\Phi}\left[\delta u_{\phi}^{(1)}\delta u_{r,t}^{(1)}+\delta u_{r}^{(1)}\delta u_{\phi,t}^{(1)}\right],$	(71)
$\displaystyle{\mathfrak{s}}_{2}({\cal P},{\cal A})=$	$\displaystyle\frac{1}{2}e^{-\Phi}\left[\delta u_{\phi}^{(1)}\delta u_{\theta,t}^{(1)}+\delta u_{\theta}^{(1)}\delta u_{\phi,t}^{(1)}\right],$	(72)
$\displaystyle{\mathfrak{S}}_{0}({\cal P},{\cal A})=$	$\displaystyle e^{\Phi}\left[{\mathfrak{s}}_{0}({\cal P},{\cal A})-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,t}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,t}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,t}\right],$	(73)
$\displaystyle{\mathfrak{S}}_{1}({\cal P},{\cal A})=$	$\displaystyle e^{\Phi}\left[{\mathfrak{s}}_{1}({\cal P},{\cal A})+\frac{2}{3}\delta S_{r\phi}^{(1)}\delta u^{\alpha(1)}_{\ ;\alpha}-\delta u^{r(1)}\delta S_{r\phi,r}^{(1)}-\delta u^{\theta(1)}\delta S_{r\phi,\theta}^{(1)}-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,r}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,r}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,r}\right],$	(74)
$\displaystyle{\mathfrak{S}}_{2}({\cal P},{\cal A})=$	$\displaystyle e^{\Phi}\left[{\mathfrak{s}}_{2}({\cal P},{\cal A})+\frac{2}{3}\delta S_{\theta\phi}^{(1)}\delta u^{\alpha(1)}_{\ ;\alpha}-\delta u^{r(1)}\delta S_{\theta\phi,r}^{(1)}-\delta u^{\theta(1)}\delta S_{\theta\phi,\theta}^{(1)}-\delta S_{r\phi}^{(1)}\delta u^{r(1)}_{\ ,\theta}-\delta S_{\theta\phi}^{(1)}\delta u^{\theta(1)}_{\ ,\theta}-\delta S_{\phi\phi}^{(1)}\delta u^{\phi(1)}_{\ ,\theta}\right],$	(75)
$\displaystyle{\cal S}_{F}({\cal P},{\cal A})=$	$\displaystyle\frac{(\varepsilon+p)e^{-2\Phi}}{\sin\theta}\bigg{\{}\frac{\delta\varepsilon^{(1)}+\delta p^{(1)}}{\varepsilon+p}\partial_{tt}{\cal Y}^{(1)}+r(\partial_{t}W^{(1)})(\partial_{tr}{\cal Y}^{(1)})+(\partial_{t}V^{(1)})(\partial_{t\theta}{\cal Y}^{(1)})+\left(r\partial_{tr}W^{(1)}+\partial_{t\theta}V^{(1)}\right)(\partial_{t}{\cal Y}^{(1)})$
	$\displaystyle+\frac{\partial_{t}(\delta\varepsilon^{(1)}+\delta p^{(1)})}{\varepsilon+p}\partial_{t}{\cal Y}^{(1)}+\left[\frac{\varepsilon^{\prime}+p^{\prime}}{\varepsilon+p}-\Phi_{,r}+\Lambda_{,r}+\frac{5}{r}\right]r(\partial_{t}W^{(1)})(\partial_{t}{\cal Y}^{(1)})+\frac{2\cos\theta}{\sin\theta}(\partial_{t}V^{(1)})(\partial_{t}{\cal Y}^{(1)})\bigg{\}},$	(76)
$\displaystyle{\cal S}_{S}({\cal P},{\cal A})=$	$\displaystyle\frac{2\mu}{r^{2}\sin^{2}\theta}\bigg{[}e^{-2\Phi}{\mathfrak{S}}_{0}({\cal P},{\cal A})-e^{-2\Lambda}\left(\partial_{r}{\cal S}_{1}({\cal P},{\cal A})+\Phi_{,r}^{\prime}\int_{0}^{t}{\cal S}_{0}({\cal P},{\cal A})dt+\Phi_{,r}{\cal S}_{0}({\cal P},{\cal A})\right)$
	$\displaystyle-\frac{1}{r^{2}}\partial_{\theta}{\cal S}_{2}({\cal P},{\cal A})-e^{-2\Lambda}\left(\frac{\mu_{,r}}{\mu}+\Phi_{,r}-\Lambda_{,r}+\frac{2}{r}\right)\left({\cal S}_{1}({\cal P},{\cal A})+\Phi_{,r}\int_{0}^{t}{\cal S}_{0}({\cal P},{\cal A})dt\right)-\frac{\cos\theta}{r^{2}\sin\theta}{\cal S}_{2}({\cal P},{\cal A})\bigg{]},$	(77)

where we note that $\delta u^{r(1)}$ and $\delta u^{\theta(1)}$ , which correspond to $W^{(1)}$ and $V^{(1)}$ , are of polar parity, while $\delta u_{\phi}^{(1)}$ , which corresponds to ${\cal Y}^{(1)}$ , is axial.

Appendix B Numerical convergence tests

In this Appendix, we examine the dependence of numerical results on the employed resolution. For this purpose, we study the time evolution of Eq. (61), assuming that ${\cal S}_{F}({\cal P},{\cal A})={\cal S}_{S}({\cal P},{\cal A})=0$ , which is equivalent to the linear perturbation analysis for the torsional oscillations. So, we should confirm that the specific oscillation frequencies obtained from fast Fourier transform (FFTs) of the time-dependent simulation data are close to the frequencies determined via the eigenvalue problem in the no-coupling limit. First, we consider a simulation with initial data given by (say)

\tilde{\cal Y}^{(2)}(r,\theta)=\frac{1}{3}\left(\frac{r}{R_{e}}\right)^{2}\left[1-\left(\frac{r}{R_{e}}\right)+\left(\frac{r}{R_{e}}\right)^{2}\right]\cos\theta,

(78)

which corresponds to an $\ell=2$ -like profile, because $\tilde{\cal Y}$ in the linear analysis is expressed as $\tilde{\cal Y}\sim\partial_{\theta}P_{\ell}/\sin\theta$ and $\partial_{\theta}P_{\ell=2}=3\sin\theta\cos\theta$ , even though the radial profile is arbitrary with the radius of the interface between the crust and envelope, $R_{e}$ .

Using the initial condition given by Eq. (78), one can expect that only $\ell=2$ torsional oscillations would be excited because the linear perturbations are not coupled with each other (recalling we have set the sources to zero). In Fig. 7, we show the FFT, using the simulation data following evolution for 1 second, by adjusting the radial grid number $N_{r}$ while keeping the angular one fixed ( $N_{\theta}=60$ ). In the same figure, we also plot, for reference, the fundamental and 1st overtone frequencies determined via the eigenvalue problem, computed as 23.9 and 898.5 Hz respectively, with dotted vertical lines. We see that the frequency of the 1st overtone determined with the simulation is severely impacted by a lower radial resolution until $N_{r}\gtrsim 500$ . To more closely compare the accuracy of the frequencies obtained from the simulation lasting one second, $f_{\rm FFT}$ (symbols), with the eigenvalues, $f_{\rm eigen}$ (dashed lines), the two values are shown as a function of $N_{r}$ in Fig. 8; the top and middle panels respectively correspond to the frequencies of the fundamental and 1st overtone torsional oscillations. Since the simulation spans 1 second, our ability to resolve individual frequencies is limited to a band of width $\sim 1$ Hz – for our purposes, such a Nyquist frequency is already sufficiently small, though we increase the run duration to 10 seconds in the next section. Nevertheless, one observes accuracy improvements with $N_{r}$ until $N_{r}\sim 500$ . To clarify this point, we also check the relative deviation, $\Delta$ , of $f_{\rm FFT}$ from $f_{\rm eigen}$ , calculated through

\Delta=\frac{|f_{\rm FFT}-f_{\rm eigen}|}{f_{\rm eigen}}.

(79)

The results are shown in the bottom panel of Fig. 8. In light of this test, we hereafter adopt $N_{r}=500$ in this study.

In a similar way, if one selects an $\ell=3$ mode-like initial condition, i.e., $\tilde{\cal Y}^{(2)}\sim 5\cos^{2}\theta-1$ , and checks the FFT using the simulation data, one can see that only ${}_{3}t_{0}$ is excited among the fundamental oscillations, while ${}_{2}t_{0}$ is not excited. Or, if one selects the sum of $\ell=2$ and $\ell=3$ mode-like profiles as an initial condition, one can see that ${}_{2}t_{0}$ and ${}_{3}t_{0}$ are simultaneously excited among the fundamental oscillations.

Next, we investigate the dependence on $N_{\theta}$ . The oscillations with larger $\ell$ have more nodes in the angular direction. So, to see the dependence on $N_{\theta}$ , it is necessary to study oscillations with large $\ell$ . For this purpose, we adopt the initial condition given by

\tilde{\cal Y}^{(2)}(r,\theta)=\left(\frac{r}{R_{e}}\right)^{2}\left[1-\left(\frac{r}{R_{e}}\right)+\left(\frac{r}{R_{e}}\right)^{2}\right]\cos^{9}\theta(1+\cos\theta),

(80)

which, importantly, does not correspond to the angular profile of any specific mode. Even so, since $\cos^{9}\theta$ and $\cos^{10}\theta$ resemble that of $\ell=10$ and 11 modes, we can expect excitations of modes corresponding to $\ell\leq 11$ . Using the initial condition given by Eq. (80) and evolving for 1 second, we compute the FFT shown in Fig. 9 for various angular resolutions: $N_{\theta}=30$ (solid), 60 (dashed), and 120 (dotted). One can observe that the deviation from the eigenvalue-determined frequencies grows with $\ell$ and is also larger for lower resolution. In Fig. 10, the relative deviation calculated with Eq. (79) is shown as a function of $N_{\theta}$ for $\ell=6-10$ . In any case, we find that one reasonably resolve modes reaching $\ell\simeq 10$ with $N_{\theta}=90$ or 120.

References

[1] N.-U. F. Bastian, D. Blaschke, T. Fischer, and G. Röpke, Universe 4, 67 (2018).
[2] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys. Rev. Lett. 50, 2066 (1983).
[3] M. Hashimoto, H. Seki, and M. Yamada, Prog. Theor. Phys. 71, 320 (1984).
[4] N. Chamel, Nuclear Physics A 773, 263 (2006).
[5] N. Andersson, Universe 7, 17 (2021).
[6] B. L. Schumaker and K. S. Thorne, MNRAS 203, 457 (1983).
[7] N. Chamel and P. Haensel, Living Rev. Relativity 11, 10 (2008).
[8] P. N. McDermott, H. M. van Horn, and C. J. Hansen, Astrophys. J. 325, 725 (1988).
[9] M. Vavoulidis, K. D. Kokkotas, and A. Stavridis, MNRAS 384, 1711 (2008).
[10] M. Ruderman, Astrophys. J. 382, 587 (1991).
[11] S. K. Lander, N. Andersson, D. Antonopoulou, and A. L. Watts, MNRAS 449, 2047 (2015).
[12] A. G. Suvorov, MNRAS 523, 4089 (2023).
[13] R. C. Duncan, Astrophys. J. Lett. 498, L45 (1998).
[14] A. L. Piro, Astrophys. J. Lett. 634, L153 (2005).
[15] Y. Levin, MNRAS 368, L35 (2006).
[16] Y. Levin, MNRAS 377, 159 (2007).
[17] C. Barat, R. I. Hayles, K. Hurley, M. Niel, et al., A&A 126, 400 (1983).
[18] T. E. Strohmayer and A.L. Watts, Astrophys. J. Lett. 632, L111 (2005).
[19] G. L. Israel, T. Belloni, L. Stella, Y. Rephaeli, et al., Astrophys. J. Lett. 628, L53 (2005).
[20] T. E. Strohmayer and A. L. Watts, Astrophys. J. 653, 593 (2006).
[21] V. Hambaryan, R. Neuhäuser, and K. D. Kokkotas, A&A 528, A45 (2011).
[22] D. Huppenkothen, A. L. Watts, and Y. Levin, Astrophys. J. 793, 129 (2014).
[23] M. C. Miller, C. Chirenti, and T. E. Strohmayer, Astrophys. J. 871, 95 (2019).
[24] H. Sotani, K. D. Kokkotas, and N. Stergioulas, MNRAS 375, 261 (2007).
[25] A. Colaiuda and K. D. Kokkotas, MNRAS 414, 3014 (2011).
[26] A. Colaiuda and K. D. Kokkotas, MNRAS 423, 811 (2012).
[27] H. Sotani, K. Nakazato, K. Iida, and K. Oyamatsu, Phys. Rev. Lett. 108, 201101 (2012).
[28] M. Gabler, Pablo Cerdá-Durán, N. Stergioulas, J. A. Font, and E. Müller, MNRAS 460, 4242 (2016).
[29] H. Sotani, K. Iida, and K. Oyamatsu, MNRAS 489, 3022 (2019).
[30] H. Sotani, Universe 10, 231 (2024).
[31] R. C. Duncan and C. Thompson, Astrophys. J. Lett. 392, L9 (1992).
[32] C. Thompson and R. C. Duncan, MNRAS 275, 255 (1995).
[33] M. Lyutikov, MNRAS 346, 540 (2003).
[34] A. M. Beloborodov and Y. Levin, Astrophys. J. Lett. 794, L24 (2014).
[35] X. Li, Y. Levin, and A. M. Beloborodov, Astrophys. J. 833, 189 (2016).
[36] S. K. Lander, Astrophys. J. Lett. 947 L16 (2023).
[37] P. D. Lasky, B. Zink, K. D. Kokkotas, and K. Glampedakis Astrophys. J. Lett. 735 L20 (2011).
[38] B. Zink, P. D. Lasky, and K. D. Kokkotas, Phys. Rev. D 85, 024030 (2012).
[39] R. Ciolfi, S. K. Lander, G. M. Manca, and L. Rezzolla Astrophys. J. Lett. 736 L6 (2011).
[40] R. Ciolfi, and L. Rezzolla Astrophys. J. 760 1 (2012).
[41] M. Gabler, P. Cerdá-Durán, J. A. Font, E. Müller, and N. Stergioulas, MNRAS 430, 1811 (2013).
[42] M. Y. Leung, A. K. L. Yip, P. C.-K. Cheong, and T. G. F. Li , Communications Physics 5, 334 (2022).
[43] R. Essick, S. Vitale, and N. N. Weinberg, Phys. Rev. D 94, 103012 (2016).
[44] M. Gabler, P. Cerdá-Durán, N. Stergioulas, J. A. Font, and E. Müller, MNRAS 443, 1416 (2014).
[45] M. Gabler, P. Cerdá-Durán, N. Stergioulas, J. A. Font, and E. Müller, MNRAS 476, 4199 (2018).
[46] H. Sotani, K. D. Kokkotas, and N. Stergioulas, A&A 676, A65 (2023).
[47] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
[48] E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein, Astrophys. J. 272, 286 (1983).
[49] H. Sotani, K. Iida, and K. Oyamatsu, MNRAS 470, 4397 (2017).
[50] C. J. Pethick and A. Y. Potekhin, Phys. Lett. B 427, 7 (1998).
[51] T. Strohmayer, H. M. van Horn, S. Ogata, H. Iyetomi, and S. Ichimaru, Astrophys. J. 375, 679 (1991).
[52] N. Andersson and G. L. Comer, Liv. Rev. Rel. 24, 3 (2021).
[53] C. Truesdell and W. Noll, The non-linear field theories of mechanics, Springer Berlin, Heidelberg (2004).
[54] L. Becerra, A. Reisenegger, J. A. Valdivia, and M. E. Gusakov, MNRAS 511, 732 (2022).
[55] A. G. Suvorov and K. Glampedakis, Phys. Rev. D 108 084006 (2023).
[56] B. Carter and H. Quintana, Proc. R. Soc. Lond. A 331, 57 (1972).
[57] S. Yoshida and U. Lee, A&A 395, 201 (2002).
[58] H. Sotani, Phys. Rev. D 107, 123025 (2023); 109, 023030 (2024).
[59] P. Pnigouras and K. D. Kokkotas, Phys. Rev. D 92, 084018 (2015).
[60] P. Pnigouras and K. D. Kokkotas, Phys. Rev. D 94, 024053 (2016).
[61] P. Arras, E. E. Flanagan, S. M. Morsink, A. K. Schenk, S. A. Teukolsky, and I. Wasserman, Astrophys. J. 591, 1129 (2003).
[62] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys. Rev. D 70, 121501(R) (2004).
[63] S. A. Teukolsky, Rev. D 61, 087501 (2000).
[64] W. Dziembowski, Acta Astronomica 32, 147 (1982).
[65] M. Gabler, U. Sperhake, and N. Andersson, Phys. Rev. D 80, 064012 (2009).