Instabilities in neutron-star post-merger remnants

Using nonlinear, fully relativistic, simulations we investigate the dynamics and gravitational-wave signature associated with instabilities in neutron star post-merger remnants. For simplified models of the remnant, we establish the presence of instability in stars with moderate $T/|W|$, the ratio between the kinetic and the gravitational potential energies. Detailed analysis of the density oscillation pattern reveals a \emph{local} instability in the inner region of the more realistic differential rotation profile. We apply Rayleigh's inflection theorem and Fj{\o}rtoft's theorem to analyze the stability criteria concluding that this \emph{inner local} instability originates from a shear instability close to the peak of the angular velocity profile and that it later evolves into a fast-rotating $m=2$ oscillation pattern. We discuss the importance of the presence of a co-rotation point in the fluid, its connection with the shear instability, and comparisons to the Rossby Wave and Papaloizou Pringle Instabilities considered in the wider literature.


I. INTRODUCTION
The merger of two neutron stars, following the gravitational-wave driven inspiral of a compact binary system, leads to the formation of a hot, differentially rotating remnant [1]. Observations of gravitational waves from such mergers are expected to shed light on the nature of matter under extreme pressures and densities by constraining the maximum mass allowed by the (hot) matter equation of state (e.g. see the recent review [2]). This puts the nonlinear dynamics of the merger into focus. It has been established that the gravitational-wave signal has robust features, most likely associated with the fundamental f-mode oscillations of the remnant [3][4][5][6][7]. Less well understood-partly because long-term post-merger simulations are prohibitively expensive (see e.g. [8])-are issues relating to the long-term survival of the remnant. The final fate of the hot remnant depends more or less directly on the involved masses (assuming that only a small amount of matter is ejected during the merger), but the timescale on which the object settles down, or collapses to form a black hole, depends on complex issues involving both dissipative mechanisms (in particular associated with neutrino emission) and magnetic field dynamics (see [9] for a recent qualitative discussion). In essence, a better understanding of merger remnant dynamics requires progress on both computational issues-how do we track the medium to longterm evolution beyond the merger?-and the physics implementation-what are the most important aspects and how do we implement them in nonlinear simulations?
This study makes a modest contribution to the discussion by focusing on the impact of the differential rotation of the remnant. This is relevant for several reasons. First, it has been established that the rotation profile of a merger remnant is quite different from that commonly assumed in studies of differentially rotating stars [10][11][12]. Most previous work on differentially rotating stars assumes that the differential rotation profile corresponds to constant specific angular momentum, leading to a profile for the angular velocity, Ω, that falls off away from the rotation axis. In contrast, merger simulations suggest that the profile should be fairly flat close to the rotation axis, rising towards a maximum at some point in the remnant beyond which it tapers off towards the profile expected for a Keplerian disk [13][14][15][16]. It seems relevant to ask how this rotation profile impacts on the dynamics, e.g. the expected f-mode oscillations. The second relevant aspect concerns the stability of these oscillations. It is well established that the oscillations of differentially rotating stars may be (dynamically) unstable already at rather modest levels of rotation (typically expressed in terms of the ratio of kinetic energy to (the magnitude of the) gravitational potential energy, T /|W |). Previous work shows that differentially rotating neutron stars become dynamically unstable when T /|W | ≥ 0.24 [17,18], but instabilities have been observed for values as small as T /|W | ≈ 0.01 for somewhat extreme rotation profiles [19]. The question then is, should we expect such low-T /|W | instabilities to be active in merger remnants? It seems a possibility worth considering, given that we do not yet have a clear understanding of the impact of the actual differential rotation profile for mergers. An initial exploration of the issue, in the context of Newtonian gravity and linear perturbations, gives an affirmative answer to the question [20]. Meanwhile, the work we present here considers the problem in full nonlinear general relativity. This is important because the differential rotation profile evolves as the system settles down, a feature that cannot be represented perturbatively. The question is if (and if so, how) this impacts on the development of an instability. Finally, a question worth asking is whether the difference in the rotation profile introduces new features not found before.
We explore these issues by carrying out numerical simulations of rapidly rotating neutron stars described by two different rotation profiles at (relatively) modest values of T /|W |, roughly at the level expected from mergers (based on the conservation of orbital angular momentum). We do not consider more extreme cases as they would mainly be of academic interest. In section II, we summarize the initial data and numerical setup for the simulations.
We also lay out the tools utilized in the data analysis. The results are reported in section III.

A. Numerical setup
We use the open source code RNSID [21] to construct the rotating initial data, representing stationary equilibrium solutions of axisymmetric relativistic neutron stars without magnetic fields. We assume that the line element for an axisymmetric and stationary relativistic space-time has the form ds 2 = −e 2ν dt 2 + e 2α (dr 2 + r 2 dθ 2 ) + e 2β r 2 sin 2 θ(dφ − ωdt) 2 , where ν, α, β, and ω are space dependent metric functions. In generating the equilibrium models, in the barotropic case, the integrability condition requires that the specific angular momentum measured by the proper time of matter is a function of Ω [22][23][24] where Ω is the angular velocity of the matter measured from infinity, v = (Ω − ω)rsinθe β−ν is the proper velocity with respect to a zero angular momentum observer. The rotation law j(Ω) used in most of previous work is the so-called j-constant law where A is a positive constant and Ω c is the angular velocity at the centre. In the Newtonian limit, j = Ω 2 , where = r sin θ is the radial distance from the rotation axis. The rotation law from Equation (3) can be rewritten as: When A → ∞, it approaches a rigid rotation, while it becomes a j-constant rotation when A → 0. That is, in this limit the specific angular momentum is constant in space (see e.g. Eriguchi and Mueller [25]). In relativity, the specific angular momentum j(Ω) is related to the metric potentials through Equation (2), so the "j-constant" law becomes: where R e is the coordinate equatorial stellar radius and the coefficientÂ = A/R e is a measure of the degree of differential rotation.
In order to investigate the dynamics of hypermassive neutron stars formed after binary neutron star mergers, the RNSID code has been modified to generate representative initial data based on a different differential rotation law introduced by Uryū et al. [26]: where p, A and B are parameters that adjust the rotation profile and j is given in equation (2) (for examples of different rotation profiles obtained from this prescription, see Figure 1 in Passamonti and Andersson [20]).

4
The initial data are then evolved using the public domain Einstein Toolkit code [27]. To evolve the fluid, we use the GRHydro module [28], together with the PPM reconstruction method and the Marquina flux formula [29]. We evolve Einstein's equations in the CCZ4 formulation [30]. A fourth-order, conservative Runge-Kutta scheme is used for the time evolution. Both the Einstein and the hydrodynamics equations are solved on a Cartesian grid using the adaptive mesh-refinement approach provided by the Carpet driver [31]. Seven levels of refinement are used to cover the simulation domain. The outermost boundary of the domain is set at ∼ 307M (≈ 460 km) and the finest refinement level has a resolution dx = 0.05M (≈ 74 m). The density of the surrounding medium (the atmosphere) has been set, relative to the initial central density ρ 0,c , to a low value ρ atm /ρ 0,c ∼ 10 −9 , leading to ρ atm ∼ 6 × 10 5 g/cm 3 . We adopt units such that c = M = G = 1 for the simulations.
Finally, we employ refluxing techniques to correct the numerical fluxes across different levels of mesh refinement [32]. The combination of the refluxing algorithm with a low atmosphere density and the CCZ4 formulation reduces position drift of the rotating profile.
This is important, as it is known that differentially rotating stars may develop spiral instabilities [33][34][35][36]. As the presence of such modes may be obscured by the numerical code not perfectly conserving linear momentum, it is important to suppress any position drift during the simulation.

B. Initial data
As our main interest is in qualitative differences and how these manifest themselves in the evolution of the system, we carry out simulations for two models. The first uses the standard j-constant rotation law with the dimensionless parameterÂ set to 1. The second model represents the rotation law (7) with p = 1,Â = 1 andB = 0.5, whereB = B/R e . In the following, we refer to this as the Uryū model to distinguish it from the j-constant case.
With this particular parameter choice we can describe the main features of the rotation laws found in merger simulations. The focus on two specific models may seem overly restrictive, but the detailed analysis we provide would not be possible for the wider parameter space.
Once we have established the tools one may consider a more exhaustive parameter survey.
We leave this for future work.
For the matter, we assume a simple polytropic equation of state, p(ρ) = Kρ Γ with Γ = 2 and K = 100. To study the effects of the rotation profile we set up two initial models with very similar bulk properties (see Table I for the main properties of these models), setting the ratio of kinetic to gravitational potential energy to a moderate value, T /|W | 0.16.
The main difference is in the rotation profile (see Figure 1). The Uryū model features a bellshaped angular velocity distribution, where the peak of the angular velocity is located at a radius of about 3.5 km. For the j-constant model the maximum angular velocity is located on the rotation axis. Note that neither of these initial models truly represents a merger remnant, as the matter distribution is truncated at a finite radius whereas a merger tends to lead to an extended disk. This should not have much impact on local features observed in the high-density region, but may affect the evolution of global dynamics. This can be tested by future work, applying our analysis of the dynamics to actual merger evolution.
in terms of the conserved density √ γu 0 ρ, where γ is the determinant of the three-metric γ ij and u µ is the fluid 4-velocity. Using the three components of the quadrupole moment in the x − y plane, we then calculate the distortion parameters η + [18,19,[37][38][39], defined as This measure serves as a proxy for the amplitude of global oscillation modes. To describe the development and saturation of the instability, we also compute the volume-integrated azimuthal density mode decomposition [34,36,40,41].
For our simulations, the numerical domain extends to ∼ 300M where the extraction of the gravitational wave signal is plausible. We extract this signal using the Newman-Penrose Scalar ψ 4 [42]. This quantity, calculated by the Einstein Toolkit module WeylScal4, is decomposed in spin-weighted spherical harmonics of spin-weight s = −2 by the Multipole module. The output is the decomposition coefficient ψ lm 4 , defined as: where dΩ stands for the differential solid angle, s Y * lm the complex conjugate of the spinweighted spherical harmonics s Y lm . The gravitational-wave strain h is linked to ψ 4 by To get the strain, we use the fixed-frequency integration method of [43] (see also Bishop and Rezzolla [44]). We choose a cutoff-frequency of 0.01 in code units, which is smaller than the initial instantaneous frequency of the waves (see Table I). To reduce boundary effects, we taper the signal using a Planck window function, with tapering width based on a period corresponding to the cutoff frequency. We also cut off the tapered parts in the final result, following the implementation in the PyCactusET module. To facilitate a comparison with detector sensitivity curves, we calculate the square root of the power spectrum density (PSD) [45]: where S h (f ) is the signal PSD, andh has been calculated as the effective strain for simplicity.

III. RESULTS AND DISCUSSION
The stability of axi-stationary fluids has been investigated in detail in the past. Whilst most instability criteria are not relevant for post-merger remnants, one possibility is the existence of a co-rotation point. Loosely this is where the pattern speed of a wave matches the angular velocity. For the particular profiles chosen here, evidence for the presence of co-rotation point instabilities was provided by [20].
All the simulations we discuss represent a star with at least one co-rotation point. We associate the instabilities that we see for these specific sets of initial data with the existence instability. As we are simulating a star and not a disk the detailed theoretical description of the precise type of instability will not carry over to our case. We have the key pre-requisites: at least one co-rotation point, with the reflecting boundary at the inner edge of the disk replaced by the centre of the neutron star, and the reflecting boundary at the outer edge of the disk replaced by the surface. However, the geometry is significantly different. In addition, the PPI, in particular, and disk physics, in general, is usually discussed at or near a j-constant state. For the Uryū case we are (at least initially) far from having constant specific angular momentum, hence the possibility of the existence of two co-rotation points.
As such, we will present analogies to the PPI and RWI cases, but will not be able to identify the instability we observe with one rather than the other.
A detailed discussion of the stability of the system will be given in section III F. In order to understand the global behaviour of the instability, we track the nonaxisymmetric modes of matter by extracting the Fourier amplitude of the density variations for the first four azimuthal multipoles, m = 1 − 4, see Equation (10). The time-dependent behaviour of the magnitude of C m illustrates the growth rate of individual modes and the nonlinear coupling associated with an instability [18], and is shown for both the Uryū model and the j-constant model in Figure 2.
As the initial data for the rotating models is close to an unperturbed axisymmetric equilibrium, at the beginning of the simulation the magnitude of each mode with azimuthal number m = 1 − 3 is essentially zero. However, the Cartesian grid naturally introduces an m = 4 mode from the numerical discretization error. The dynamics of this m = 4 mode is similar to previous findings in the literature [18,34,39]. For this reason we consider the m = 4 mode to be "grid noise" in these simulations. Dynamically unstable rotating stars, as they deform to non-axisymmetric configurations, may be relevant sources of quasi-periodic gravitational waves (see e.g. Shibata et al. [38]).
The gravitational-wave signature depends on the specific evolution of the leading quadrupole mode and the extent to which other modes suppress its growth. Baiotti et al. [18] find that generic nonlinear mode-coupling effects appear during the development of the instability and these can severely limit the persistence of any bar-mode deformation. Despite the initial data used here having a T /|W | well below that needed for the classical bar mode instability, we note the growth of low-m modes in a similar fashion to studies of bar modes such as Baiotti et al. [18]. Figure 3 shows the gravitational-wave strain as viewed by an observer located on the z-axis at a distance of 100 Mpc. This figure should be viewed alongside Figure   2, as there is an association between gravitational-wave pattern and oscillation modes, which can also be found in previous work (e.g. [39]).
. Both profiles show a dominant peak at nearly the same frequency. For the Uryū model, the peak is at 2.9 kHz, while the peak is located at at 3.0 kHz for the j-constant model.
The main features are the same in both models. As the instability grows, a long-lived gravitational wave of modest amplitude is generated in both polarization. As shown in the previous section, the instability is dominated by an m = 2 mode throughout the simulation, which in turn dominates the gravitational wave emission. We distinguish some modulation of the signal on long timescales in the j-constant model. By comparing the gravitational wave emission in Figure 3b to the detailed mode growth in Figure 2b it could be argued that this modulation is linked to the m = 2 mode behaviour. However, as all modes appear to have saturated by the end of simulation for both models, this suggests that the gravitational wave emission will remain largely unchanged for both models, and so the two models will be difficult to distinguish.
The power spectra associated with the gravitational-wave signals are shown in Figure 4, comparing the Uryū (left) and j-contant (right) models. The main peak, located at around 3 kHz, is very similar for the two rotation profiles. We now consider the dynamics associated with the mode instability, comparing and contrasting the behaviour for the Uryū and j-constant models. In Figures 5 and 6 the snapshots of the normalized density deviation from the averaged value, (ρ − ρ )/ ρ , where . . . represents an angular average at a specific radius in the equatorial plane [46].
The computational domain is divided into two parts, approximately representing the neutron star interior and the surrounding low-density medium (the atmosphere). Separate colour maps are used to visualize the density variation inside, and outside of, the neutron star.
In the case of the Uryū rotation law, the initial density variation originates from the grid discretization. This leads to a quasi-stationary m = 4 mode pattern, evident in the first (t = 1 ms) snapshot in

D. Rotation profile evolution and co-rotation radius
Given a specific oscillation mode, we can associate the frequency to the angular velocity of the rotating profile, and then define the co-rotation point as the position where the mode's pattern speed matches the bulk angular velocity. It has been suggested (see, e.g., [20,39,47]) that a low T /|W | instability sets in when such co-rotation points of the unstable modes exist. For the j-constant rotation law, only one co-rotation point can possibly exist for each mode. When we consider the Uryū model, however, a given mode may exhibit two distinct co-rotation points due to the bell-shaped feature of the angular velocity profile (see Figure 1 and the discussion in [20]).
In Figure 7, we present the time evolution of the angular velocity profile for both models. with the m = 2 mode locates at around 4 km, and does not change with time.
The additional co-rotation point within the Uryū model makes the mode analysis more complicated, leading us to consider the relevance of each co-rotation point. Will there, for example, be an independent mode pattern associated with the inner co-rotation point?
Based on the results of the density deviation plots for the Uryū model, we find that the oscillation pattern close to the center rotates faster throughout the simulation. To analyze this inner local oscillation, we calculate the distortion parameter η + (see Equation (9) Figure 8. Note that the inner local oscillation, inferred from Figure 5, is located within a radius of 5 km. The first region then represents the inner core. The second region bridges the core with the outer envelope. The third and fourth regions represent the envelope of the remnant. We calculate the Fast Fourier Transform (FFT) of the time series of the distortion parameter η + inside each region. We divide the time series into windows with a width of 4.5 ms to perform the FFT calculation.
The power spectra are shown in Figure 9.

E. The inner local instability
In this section, we discuss the nature of the inner fast-rotating local oscillation for the Uryū model. To minimize the grid effects, we conduct two further simulations for both models without the innermost refinement level, which has the boundary set at around 3.5 km (close to the peak of the angular velocity profile in the Uryū case). The innermost refinement boundary for these new simulations extends to a radius of about 6 km, covering the entire peak region of the angular velocity for the Uryū model. Reassuringly, this does not lead to qualitative changes to the results, but it nevertheless ensures that the grid has no impact on the inner fast-rotating local oscillation.
We then calculate the xy-component of the vorticity 2-form [40], in the equatorial plane. As a reminder, h = 1 + + P/ρ 0 represents the specific enthalpy, where is the internal specific energy and P the pressure. The time evolution plot of the Ω xy contour for the Uryū and j-constant models is shown in Figures 10 and 11, respectively.  At this point the whole core region deforms into a bar-shape. For the j-constant model (see Figure 11), a persistent spiral pattern develops around the centre throughout the simulation.
These different instability behaviours distinguish the Uryū model from the j-constant case.
It has been found that differentially rotating barotropic flows with sufficiently strong shear layers are unstable to the formation of vortex chains (see [48][49][50][51] for examples). The Uryū model features a ring-shaped flow that rotates relatively faster than the rest of the remnant, in analogy to a jet flow within a rotating fluid (see e.g. Solomon et al. [51]). We suspect a similar shear instability mechanism plays a role in the observed formation of vortices. In the following, we will test an instability theorem to get a glimpse of the underlying physics behind the inner local instability.

19
FIG. 11: Same as Figure 10, but for the j-constant model. The dotted circle is located at the same radius as in Figure 10. Compared to the Uryū model, a persistent rotating spiral pattern develops around the center.

F. The instability criteria for rotating flow
In the study of flow instabilities, Rayleigh first developed a general stability theory for inviscid parallel shear flows, and showed that a necessary condition for linear instability is that the velocity profile has a point of inflection [52]. Fjørtoft then gave a strict necessary condition that there is a maximum of vorticity for inviscid instability [53]. For the stability of inviscid rotating flows, Rayleigh also obtained a criterion which is the analogue of the inflection point theorem in parallel flow [54].
The detailed derivation of the generalization of Fjørtoft's theorem to rotating flow can be found in [55]. It states that a necessary, but not sufficient, condition for linear instability of inviscid rotating flow is that ξ (Ω − Ω s ) < 0 somewhere in the flow field. Here is the vorticity of the background flow, Ω(r) is the mean angular velocity, and Ω s = Ω(r s ) is  Figure 12 shows the instability criteria results in the equatorial plane. We find the instability criteria also extend to the vertical plane. We plot the angular velocity contour in the − z plane, and indicate the locations for the inflection points and the maximum angular velocity along the direction (as in Figure 13). As we can see, at each horizontal plane with different values of z, there exists an inflection point in the region inside the maximum angular velocity. Instabilities similar to what we see in Figure 10 should thus be expected to occur at each level of z. In analogy to Taylor-Couette flow, the instability is essentially three dimensional. In Figures 14 -15, we plot the time snapshots of vorticity Ω yz in the y-z plane for the Uryū and j-constant models.
In Figure 14, disrupted. We suspect that by this time the instability reaches the nonlinear regime (which is also suggested by Figure 10). Figure 15 shows the Ω yz plot for the j-constant model.
In comparison with the Uryū model, we do not see clear evidence for a linear instability throughout the early phase of the time evolution in the Ω yz plot. This is in agreement with our previous stability criteria results which suggest that the flow in the j-constant model is stable.
G. Low-T /|W | instability and co-rotation points Now let us switch our attention to the general low-T /|W | instability. The linear analysis of Watts et al. [47] suggested that low-T /|W | instabilities are triggered when the co-rotating f -mode enters the co-rotation band within the differentially rotating star. By investigating the distribution of the canonical angular momentum, Saijo and Yoshida [56] found that the instability sets in around the co-rotation radius of the star, and grows as there is an inflow of angular momentum inside the co-rotation radius. In this picture, the co-rotation radius has a crucial role in that a wave propagating radially across it can be amplified. However, significant growth of a wave typically requires many passages through the co-rotation point [57]. It appears as though a resonant cavity is required to drive the modes in co-rotation to large amplitude [33]. For example, in the case of the Papaloizou-Pringle instability (PPI), the inner and outer edges of the disk or torus forms a resonant cavity in which waves are reflected back and forth [58].
Similarly, Lovelace et al. [59] have analyzed the so-called Rossby wave instability (RWI) in Keplerian accretion disks and found that it occurs when there is an extremum in the radial profile of L(r) ≡ (ΣΩ/κ 2 )S 2/Γ , where Σ is the surface mass density of the disk, Ω is the angular rotation rate, S(r) is the specific entropy, Γ is the adiabatic index, and κ is the radial epicyclic frequency. Extrema of L(r) could come from several sources. In [59] they considered the special case where there is a local maximum in the disk entropy profile, S(r).
This maximum acts to trap the waves in the vicinity of the maximum, given sufficiently strong variation. Li et al. [57] have presented a detailed linear theory for the RWI and show that it exists for a wider range of conditions, specifically, for the case where there is a "jump" over some range of r in Σ(r) or in the pressure P (r). They also point out that the profiles of Σ(r) and P (r) considered are not the only ones which may lead to instability. For example, a profile with local extreme in the vortensity distribution (see Equation (18)

25
Closely related to this, in the study of differentially rotating neutron stars, Ou and Tohline [33] show that models of differentially rotating neutron stars can also exhibit a local minimum in their radial vortensity profile, and a similar resonant cavity mechanism seems to trigger the one-armed spiral instability. In addition to the one-armed (m = 1) spiral mode, they have found that higher order (m = 2 and 3) nonaxisymmetric modes can also become unstable if the associated co-rotation points that resonate with the eigenfrequencies of these higher modes also appear inside the star. Note that their model configurations feature centrally condensed, rather than toroidal, density structures. These studies suggested that the presence of a minimum in the profile of the vortensity (see Equation (16)) is a necessary condition for a mode in co-rotation to be unstable [33,39].
Given the different shape of the angular velocity for the Uryū model and j-constant model, we expect the vortensity profile for the two models to be different. To verify this and facilitate further analysis, we have computed for both models the Newtonian vortensity, defined as the ratio, along the radial cylindrical coordinate, between the radial vorticity and the density [39,57], i.e.
where κ 2 = 2Ω r d dr (r 2 Ω) is the square of the radial epicyclic frequency (so that κ 2 /2Ω is the vorticity) [62]. The results are shown in Figure 16. The co-rotation points of the global m = 2 mode (see Figure 7) are also indicated. Both co-rotation points are located near the minimum of each vortensity profile. This is in agreement with previous studies (see e.g. [39]). Now let us focus on the inner local oscillation for the Uryū model. Based on Figure 9, the oscillation has an frequency of 3.6 kHz and the dominant mode azimuthal number is m = 2 (inferred from Figure 5) during the time period t > 10ms. With a pattern speed of 1.8 kHz, its co-rotation point is close to the peak of the angular velocity profile (see Figure 7).
During the period where the oscillation develops (roughly t = 8 ms − 16 ms), the angular velocity profile of the remnant changes. Still, we can roughly find the "co-rotation region" of the local mode which is close to the peak of the vortensity profile (see Figure 16).

26
H. Linear and nonlinear instability for the inner local instability Figure 10 illustrates the process of the formation and development of the inner local instability within the star's core. Roughly speaking, it goes through three stages. First, the rapid growth of the initial small amplitude perturbations. Next follows the production of large-scale vortices and their interactions with the background flow. Finally there is a coupling of vortices with global spiral waves. For the first stage, roughly t < 8 ms, the shear instability described in subsection III F plays a major role. Later on, a nonlinear disturbance leads to the formation of vortices and the oscillation pattern. These stages are similar to the case of the RWI in thin accretion disks with density or pressure structures (see e.g. [63]). When we discuss the inner local oscillation for the Uryū model, the shear instability based on Rayleigh's and Fjørtoft theorem, and co-rotation resonance effects from RWI/PPI may be well entangled. At the end of the simulation, the inner local oscillation almost synchronises with the outer f -mode oscillation.

IV. CONCLUSIONS
We have carried out numerical simulations of rapidly and differentially rotating neutron star configurations, inspired by post-merger remnants. The results demonstrate that different angular velocity profiles lead to slightly different dynamics, impacting on the growth and saturation of the low-T /|W | instability. In particular, we find that the profile of the Uryū model generates a more dominant m = 2 mode perturbation. For the j-constant model with similar bulk properties, the m = 2 mode is more strongly coupled to other multipole modes (especially an m = 3 component). In this case, the mode coupling generates distinct gravitational-wave bursts rather than a continuously growing amplitude in the beginning.
In an addition, we find the oscillation pattern close to the centre of the remnant behaves differently in the Uryū model, a feature that triggers a local instability. We show that this local instability is directly linked to the bell-shaped feature of the angular velocity profile, and occurs in the inner part of the remnant where strong shear layer exists. We apply the generalized Fjørtoft's theorem to the rotating profiles, and find that the Uryū model has inflection points in its angular velocity profile which satisfy the instability criteria, while the j-constant model appears to be stable according to this measure. The vorticity contour of the Uryū model confirms that the inner local instability occurs in the predicted unstable region. This inner local instability starts with a linear instability in the horizontal as well as in the vertical plane. It then leads to the formation of vortices, which merge together to form a fast rotating m=2 oscillation, distinguishable from the global m=2 mode that appears in the rest of the remnant. As time goes on, the inner local m = 2 oscillation synchronizes with the global m = 2 mode. For the j-constant stellar model, we only observe the global mode development.
For the general low-T /|W | instability, we find that an m = 2 f-mode co-rotation point exists inside the rotating profile in both cases. This co-rotation point is located near the minimum of the corresponding vortensity profile. This indicates that a co-rotation resonance may amplify the magnitude of the f-mode as discussed in the case of RWI/PPI. For the Uryū model, the co-rotation point for the local m = 2 oscillation is located near the peak of the angular velocity profile. The linear shear instability and RWI/PPI may well participate in the development of this inner local oscillation.
This study provides an initial survey of the nonlinear effects associated with unstable modes for different rotation laws, complementing the linear perturbation study from [20].
We focused on comparing a rotation profile inspired in binary neutron star merger remnants to the standard j-constant rotation law. The results provide qualitative insights into the impact of the rotation profile on the development of mode instabilities. We have performed some additional simulations with the same models, but with additional parameters chosen so that no co-rotation points exist, and these show no instability. However, we have not found the instability threshold to high accuracy. A more detailed parameter survey, exploring the dependence on the different parameters, like the peak and position of the angular velocity, may lead to a deeper understanding about the instability and mode dynamics during the post-merger phase. The discovery of the inner local instability highlights the importance of the study of instabilities of shear flows in the framework of relativity. Also, the effects of other important ingredients like realistic (hot) equations of state, magnetic fields and the presence of matter in a disk surrounding the remnant on the dynamics and the instabilities require further investigation.
The first release of the initial data, parameter file, analysis and visualization scripts is available through Zenodo [64]. 28