Continuous gravitational waves from neutron stars: current status and prospects

Gravitational waves astronomy allows us to study objects and events invisible in electromagnetic waves. It is crucial to validate the theories and models of the most mysterious and extreme matter in the Universe: the neutron stars. In addition to inspirals and mergers of neutrons stars, there are currently a few proposed mechanisms that can trigger radiation of long-lasting gravitational radiation from neutron stars, such as e.g., elastically and/or magnetically driven deformations: mountains on the stellar surface supported by the elastic strain or magnetic field, free precession, or unstable oscillation modes (e.g., the r-modes). The astrophysical motivation for continuous gravitational waves searches, current LIGO and Virgo strategies of data analysis and prospects are reviewed in this work.


Introduction
Gravitational-wave (GW) astronomy has been one of the fastest-growing fields in astrophysics since the first historical detection of a binary black-hole (BH) system GW150814 (Abbott et al., 2016a). In addition to studying the nature of gravitation itself, it may be used to infer information about the astrophysical sources emitting the GWs. This review concentrates on a specific kind of prospective GWs: persistent (continuous) gravitational waves (CGWs), emitted by neutron stars (NSs). The article is arranged as follows. Section 1 gathers introductory material: Section 1.1 presents the basics of the GWs theory, Section 1.2 contains a brief overview of GWs detections, Section 1.3 describes properties of NSs and features of hitherto detected NSs-related GWs-a binary NS merger GW170817 (Abbott et al., 2017d), Section 1.4 gathers general information about CGWs, whereas Section 1.5 is devoted to the main data analysis methods used in CGWs searches. Following sections describe the main CGWs emission mechanisms: elastic deformations (Section 2), magnetic field (Section 3), oscillations (Section 4), free precession (Section 5). Finally, in Section 6 contains summary and discussion.

Basics of the Gravitational Radiation Theory
According to the general theory of relativity (Einstein, 1916(Einstein, , 1918, GWs are perturbations in the curvature of space-time, travelling with the speed of light. To produce waves, just as in the case of electromagnetic (EM) waves, accelerated movement of charges (masses) is needed. The lowest allowed multipole is the quadrupole, as the monopole is forbidden by the mass conservation and the dipole by the momentum conservation. A non-negligible time-varying quadrupole moment may be provided by e.g., binary systems: BHs or NSs, rotating non-axisymmetric objects (i.e., deformed NSs) or non-spherical explosions (supernovae). According to the quadrupole formula at the lowest order (Einstein, 1916(Einstein, , 1918, GW amplitude strain tensor h i j at position r is is the mass-quadrupole moment in the transverse-traceless (TT) gauge 1 , evaluated at the retarded time (t − r/c) and ρ is the matter density in a volume element d 3 x, at the position x i ; c and G is the speed of light and gravitational constant, respectively. Extension of the above h i j expression includes second-order multipole moment, called current-quadrupole moment, given by Thorne (1980). It describes the dynamics of the mass currents that can lead to GWs emission caused by e.g., the r-mode instability (discussed in Section 4).
Equation (1) shows that an axisymmetric NS rotating along its axis will not emit GWs, because its mass-quadrupole moment will not vary in time. The GW luminosity L GW is with ... brackets denoting time averaging, and a dimensionless parameter quantifying the level of asymmetry. The moment of inertia along the rotational axis I 3 scales with NS mass M and radius R as I 3 ∼ MR 2 ; ν is NS rotational frequency. An estimate of the GW strain amplitude is thus h 0 ∼ 10 2 G I 3 ν 2 c 4 d , which is inversely proportional to the distance to source d. Propagation of the GWs in vacuum is governed by a standard wave equation: for which the simplest solution is a plane wave solution: where k is a wave 3-vector defining the propagation direction, related to the wavelength λ as λ|k| = 2π, x is a 3-vector of coordinates, A i j is constant amplitude and α (i)( j) is the constant initial phase. ω GW = 2π f GW is rotational (angular) frequency, while f GW is a frequency of the gravitational wave.
In the TT gauge, the above equation can be rewritten in the following form that depends on two independent polarisations of the GW: plus '+' and cross '×' (see explanation in Jaranowski & Królak 2009): Here h + and h × are the polarisations of a plane wave moving in the +z direction: Full derivation can be found e.g., in Misner et al. (1973); Bonazzola & Gourgoulhon (1996); Jaranowski & Królak (2009);Prix (2009).
Note that h + and h × depend on the mechanism of gravitational radiation and will be discussed later in this review, while F + and F × are periodic functions with the period equal to one sidereal day (for ground-based detectors such as the LIGO or Virgo), due to the rotation of the Earth. Additionally, F + and F × depend on the wave polarisation angle ψ, the angle between detector's arms ζ (usually ζ = π/2), and two amplitude modulation functions a(t) and b(t), which depend on the location and orientation of the detectors on Earth and the position of the GWs source on the sky (full representation of a(t) and b(t) can be found e.g., in Jaranowski & Królak 2009, Appendix C). Order of magnitude of frequency variations due to the daily and annual motion of the Earth is thus 10 −6 Hz and 10 −4 Hz, respectively.

Brief History of Gravitational Waves Detections
First indirect evidence of the existence of GWs was deduced from the observations of stars: binary systems containing white dwarfs (Paczyński, 1967), and later neutron stars, most notably the Hulse-Taylor pulsar (catalogue number: PSR B1913+16, Hulse & Taylor 1975. The loss of orbital energy manifests itself as the shrinkage of the orbit and in the result as a drop of the orbital period (Peters & Mathews, 1963): (1−e 2 ) −7/2 1 + 73 24 e 2 + 37 96 e 4 M 1 M 2 (M 1 + M 2 ) −1/3 , (11) where P-orbital period, e-orbital eccentricity, M 1 and M 2 -masses of the components. In 30 years of observations  observed orbital decay was consistent up to 0.13 ± 0.21% level with the theoretical prediction for the emission of GWs (Weisberg & Taylor, 2005). The observed orbital decay (Taylor & Weisberg, 1982, 1989Weisberg et al., 2010) is in an excellent agreement with the expected loss of energy due to the GW radiation as described by the general relativity.
In 2015, the real GW astronomy began with the first direct detection of the binary BH coalescence, GW150914 (Abbott et al., 2016a). The signal was registered by two LIGO interferometers (Aasi et al., 2015b) and analysed jointly by the LIGO and Virgo Collaborations (LVC). Thus far, nine more binary BH mergers (Abbott et al., 2016b(Abbott et al., , 2017a(Abbott et al., ,b,c, 2018b were reported. Virgo Collaboration operates a third detector of the global network (Acernese et al., 2014). A second breakthrough came with the LVC observation of the binary NS merger (Abbott et al., 2017d). A network of three GW detectors cooperated with many EM observatories, performing first observations of GWs and a broad spectrum of EM waves from the same source (Abbott et al., 2017e,f). These unique, multi-messenger observations allowed for the first 'standard siren' measurement of a Hubble constant (Abbott et al., 2017g), measurement of the GW propagation speed (Abbott et al., 2017d), discovery and study of the closest kilonova event (Abbott et al., 2017h), estimation of the progenitor properties (Abbott et al., 2017d(Abbott et al., ,i, 2019a and study of the post-merger remnant (Abbott et al., 2017j, 2019h).
Recently, LVC published the first catalogue with all the previous detections and source parameters' estimates 2 .
At the time of writing, second-generation of interferometers: the Advanced LIGO (aLIGO, Harry et al. 2010) and Advanced Virgo (AdV, Acernese et al. 2014) -are collecting data. 'Second-generation' refers to the strongly improved versions of the initial, first-generation detectors. First-generation observatories were the following GW interferometers: TAMA (300-m arms) near Tokyo, Japan (Takahashi, 2004), GEO600 (600-m arms) near Hannover, Germany (Willke, 2002), Virgo (3-km arms) near Pisa, Italy (Freise et al., 2005), and LIGO (two instruments with 4-km arms each) in Hanford and Livingston, US (Abramovici et al., 1992). Soon, next detectors will join the global GW network: the Japanese collaboration that built TAMA is now testing the 2nd-generation underground and cryogenic detector KAGRA-KAmioka GRAvitational wave telescope (Akatsu et al., 2017) and the third LIGO interferometer will be placed in India and operated by the Indian Initiative in Gravitational-wave Observations (IndIGO) in near future (Unnikrishnan, 2013). Next step will be to design, build and operate third-generation detectors, with the planned sensitivity an order of magnitude better than the second-generation detectors. A European consortium is in the design stage of the Einstein Telescope (ET)-the underground trio of triangular Figure 1: Comparison of the characteristic strains as a function of GW frequency for existing and planned ground-based detectors (for details see Moore et al. 2015; plot generated with an interactive tool at http://gwplotter.com). The expected CGWs amplitudes ranges of pulsars are marked with the pastel-blue region. For comparison, the GW150914 signal strain is represented by the pink outline.
interferometers, each 10 km long (Sathyaprakash et al., 2012). Comparison of the current and future ground-based detectors characteristic strains in the context of CGWs sensitivities (denoted as 'Pulsars' on the plot) is shown in Figure 1.
Underground interferometers will probe GW frequencies down to ∼1 Hz, but to reach even lower frequencies (astrophysically interesting from the point of view of e.g., extreme mass ratio inspirals, heavy binary systems inspiral, or primordial GWs fluctuations from the early Universe), space-based interferometry is required. A visionary project LISA (Laser Interferometer Space Antenna), led by European Space Agency includes a triangular configuration (each arm 2.5 × 10 6 km long) of three satellites that will be placed in a solar orbit at the same distance from the Sun as the Earth (Amaro-Seoane et al., 2017). The mean linear distance between the formation and the Earth will be 5 × 10 7 km.
Currently, ongoing upgrades of the existing detectors, new detectors in the network, and improvements of the data analysis methods will lead to an increase in sensitivity and better quality of the detector data and, as a result, to the detections of the weaker signals. GW sources may be divided in to four categories, depending on the duration and strength of the signal, as shown in Table ??: continuous gravitational waves (CGWs), subject to this review; stochastic background GW which is a mixture of a large number of independent sources; inspirals and mergers of binary systems; burst sources e.g., supernovae explosions or magnetar flares. Search strategies for each type  depend on the duration of GWs emission, knowledge of the signal model, characteristic amplitudes of the GWs, and the computational resources. So far only compact objects' inspirals were detected in the LIGO-Virgo observational runs. It is expected that in the next LIGO and Virgo observing seasons, more subtle signals, such as CGWs or a stochastic background, will be detected. Steady improvement of the search methods and sensitivity of the detectors was demonstrated in the past: searches for CGWs (Abbott et al., 2017k, 2018c(Abbott et al., 2017k, , 2019b, stochastic background (Abbott et al., 2017l, 2018d and burst signals (Abbott et al., 2017m, 2019d. Even though no significant detections were claimed, astrophysically interesting upper limits were determined. They will be discussed in the next sections.

Properties of Neutron Stars
Pioneering hypothesis of existence of the dense stars, that look like giant atomic nuclei, was given by Landau (1932), even before the discovery of a neutron 3 . Existence of the NSs, as remnants of the supernovae explosions, was proposed by Baade & Zwicky (1934), just after the breakthrough discovery of a neutron by Chadwick (1932). This hypothesis waited more than 30 years to be confirmed, when Jocelyn Bell Burnell first discovered a 'rapidly pulsating radio source' (Hewish et al., 1968a,b), which was interpreted subsequently as a fast-spinning NS that emits a beam of electromagnetic radiation similarly to a lighthouse. Since then, the nature of NSs is still a mystery. The extremely large densities and pressures present inside NSs cannot be reproduced and tested in terrestrial laboratories. We also do not have a credible theoretical description of how matter behaves above the saturation density of nuclear matter n 0 = 0.16 fm −3 (density at which energy per nucleon of infinite, symmetric nuclear matter has a minimum).
By measurements of the NSs masses M and radii R one can, in principle, determine properties of the matter inside NSs, such as the relation between its pressure P and density ρ (the currently unknown equation of state of dense matter, EOS; for a textbook introduction to the subject, see Haensel et al. 2007). Bijective functions P(ρ) and M(R) are schematically represented on Fig. 2.
Putting constraints on EOS from observations can be done using the fact that M(R) is a bijective function to P(ρ), typically by solving hydrostatic equilibrium for a spherically symmetric distribution of mass, named Tolman-Oppenheimer-Volkoff (TOV) equation (Tolman, 1939;Oppenheimer & Volkoff, 1939): constrained by the equation for the total gravitational mass inside the radius r: Descriptions of the various methods used in research on EOS determination can be found in e.g., the review by Özel & Freire (2016). From the EM observations we know that any realistic EOS should support NSs with masses around 2 solar masses (Demorest et al., 2010;Antoniadis et al., 2013;Fonseca et al., 2016). Uncertainties in the NSs radii and crust properties reflecting our limited knowledge of the EOS are discussed e.g., by Fortin et al. (2016).
Although this review focuses on the long-duration CGWs, for the sake of completeness, in the remainder of this section, we will provide a brief description of the tidal deformability effect imprinted on the GWs emitted during the last orbits of a binary NS coalescence. So far, one measurement of this kind-the GW170817 event-was successfully performed and published (Abbott et al., 2017d(Abbott et al., , 2018f, 2019a, whereas a second event, dubbed S190425z, was recently reported by the LVC via the public database service (LIGO/Virgo GraceDB, 2019). For the recent theoretical studies concerning the interpretation of the GW170817 tidal deformability measurement and its relation to the dense-matter EOS, see Annala et al. (2018); Burgio et al. In a binary system, a quadrupole moment of each NS is induced by the companion NS, due to the presence of the external tidal field ε i j (Misner et al., 1973;Hinderer, 2008): where the proportionality factor λ td is called the tidal deformability parameter, expressed in the lowest order approximation as The parameter k 2 is the quadrupole (l = 2) tidal Love number (Love, 1911;Flanagan & Hinderer, 2008): with the compactness of object x = GM/Rc 2 (M denoting the gravitational mass, Rradius), and y-a solution of evaluated at the stellar surface r = R (Flanagan & Hinderer, 2008;Van Oeveren & Friedman, 2017). For convenience, a dimensionless value of the tidal deformability Λ is often defined as Λ = λ td GM/c 2 −5 .
Tidal contribution to the GW signal are extracted from the last stages of the inspiral phase. What is actually measured is, in the first approximation, the effective tidal deformabilityΛ, a mass-weighted combination of individual dimensionless tidal deformabilities Λ 1 , Λ 2 , defined as with M 1 , M 2 denoting the component gravitational masses 4 . In principle, by solving the above equations, one can relate the measurable value ofΛ to the EOS parameters P and ρ.
4 Note that the most precisely measured quantity from an inspiral phase is a so-called chirp mass: M = (M 1 M 2 ) 3/5 (M 1 +M 2 ) 1/5 . Using Newton laws of motion, Newton universal law of gravitation, and Einstein quadrupole formula, one can see that M depends only on GW frequency f GW and its derivativeḟ GW -

General Information about Continuous Gravitational Waves
According to the ATNF (Australia Telescope National Facility) Pulsar Database 5 (Manchester et al., 2005), more than 2700 pulsars are known. Assuming that CGW signal primarily corresponds to twice the NS spin frequency (in case of elastic deformations, see Section 2) or close to 4/3 of the spin frequency (Newtonian approximation in case of r-modes, see Section 4), around 300 pulsars from ATNF Database are in the LIGO and Virgo detectors' frequency range which is the most sensitive around 50-800 Hz, as shown in Figure 1. Around 200 of these pulsars have precise ephemerides, as well as, measured ν andν. We also know from evolutionary arguments that the Galaxy contains ∼10 8 -10 9 NSs (Narayan, 1987;Camenzind, 2005), of which ∼160,000 are isolated (Lorimer, 2005). For recent results concerning the population synthesis of isolated radio pulsars in the Galaxy see Popov & Prokhorov (2007); Lorimer (2011);Cieślar et al. (2018).
As was observed by Patruno et al. (2017), observed limit of NSs spins is currently ν max ≈ 700 Hz. Two possibilities were considered in Haskell et al. (2018a): (a) ν max corresponds to a maximal allowed spin, above which the centrifugal forces causes mass shedding and destroy the star (also called the Keplerian frequency). As was shown in Haskell et al. (2018a), ν max cannot be less than ≈ 1200 Hz, and that the observed lack of NSs spinning faster than ≈ 700 Hz is not consistent with minimal physical assumptions on hadronic physics; (b) presence of additional spin-down torques acting on the NSs, possibly CGWs emission.
It means that in our Galaxy numerous promising CGWs sources are hidden and awaiting for detections of their gravitational signatures.
Rotating non-axisymmetric NSs (Ostriker & Gunn, 1969;Melosh, 1969;Chau, 1970;Press & Thorne, 1972;Zimmermann, 1978) are expected to emit CGWs due to the existence of time-varying mass-quadrupole moment. Such signals have smaller characteristic amplitudes than signals emitted from compact binary mergers, but their overall duration is much longer. In the case of almost-monochromatic CGWs, their integrated signal-to-noise ratio (S NR) increases with the observation time T as where h 0 is the instantaneous GW strain amplitude and S is the amplitude spectral density of the detector's data signals' frequency. For comparison, GW150914 lasted for T ∼ 0.2 s in the sensitive part of detectors' band with the average GW amplitude h 0 ∼10 −21 , yielding S NR∼24. For CGWs, the expected amplitude is a few order of magnitudes smaller, h 0 ∼10 −25 , but the data collection lasts for T of the order of months quantities determined directly from the observational data: M = c 3 G 5 96 3 π −8 f −11 or even years. This is one of the incentives for consideration of the CGWs as serious candidates for the future detections (Brady et al., 1998;Jaranowski & Królak, 2000). By simple manipulation with Equation (20), one can roughly estimate the observational time that is needed to observe signals considered in this paper. Let us put T ∝ S (S NR/h 0 ) 2 . While writing this paper, currently ongoing O3 observational run takes place and LIGO detector reaches √ S ∼ 10 −23 Hz −1/2 (see Abbott et al. 2019i also for comparison of sensitivities of Virgo and KAGRA detectors, as well as predictions for the future observing seasons). Typically, a threshold for the CGW detection is set S NR th = 5. For the most promising scenario of CGW emission (triaxial ellipsoid, see Section 2), h 0 ∼ 10 −26 (see Equation (27)). As a result, in order to detect such a signal, one needs T ∼ 300 days of good quality, coherent data (note that O3 run is scheduled for 1 year, what makes the future CGW searches very promising). Analogously for the magnetic field distortions model (Equation (35)), where h 0 ∼ 10 −30 , almost 80 mln years of data is needed!

Methods and Strategies of CGWs Searches
Detectability of the CGWs signals depends on the observational time (Equation (20)), but also on the balance between computational cost of the accurate data analysis and computational resources. For some isolated NSs the relevant parameters, such as sky position (e.g., right ascension α and declination δ), rotational frequency ν ( f GW is, in case of an elastic deformation, a mixture of 1ν and 2ν, see Section 2, or f GW = 4/3ν for r-modes, see Section 4) and spin-downν are known from EM observations 6 . For these objects, dedicated targeted searches are performed, in order to check if a CGW signal is associated with known parameters of the individual pulsar (Nieder et al., 2019). Such searches are computationally easy to perform. A slight modification to the targeted searches are so-called narrow-band searches, which allow for a small mismatch between the frequency parameters known from the EM observations and the expected GW signal.
Another type of search strategy, in which position of the source in the sky is assumed, but the frequency parameters are unknown, is called a directed search. It may be applicable to e.g., a core-collapse supernovae remnants. Example of a young and relatively close supernova remnant, for which spin frequency (and its derivatives) are unknown, is the Cassiopeia A remnant (Fesen et al., 2006). As was shown in Wette et al. (2008), CGW strain depends on the assumed age a and distance d. Additional intrinsic parameters such as NS mass or its equation of state increase uncertainty for the GW strain of the strongest possible signal h age 0 , by 50%: Additionally, young and hot NSs may become unstable and undergo various dynamical processes (e.g., cooling processes and oscillations-see Section 4). Full realistic description of the CGWs emission, in presence of multiple physical phenomena may require inclusion of higher frequency derivatives and in general is extremely hard to model. Nevertheless, searches were performed for the CGWs from known supernovae remnants (Abadie et al., 2010;Zhu et al., 2016;Ming et al., 2016Ming et al., , 2018Ming et al., , 2019Abbott et al., 2019h), as well as in the direction towards the Galactic Centre region, where massive stars (progenitors of supernovae explosions) are found in stellar clusters (Aasi et al., 2013;Dergachev & Papa, 2019c).
The most computationally intensive searches are all-sky (blind) searches, since only minimal assumptions are made to search for the signals from a priori unknown sources. Such a search requires well-optimised and robust tools, because the less is known about the source, the smaller sensitivity of the search can be achieved and bigger computational cost is required. All types of searches are summarised in Figure 3. Several methods (described below) were so far tested and used in blind searches on mock and real data. Each all-sky search has different advantages in the sensitivity vs. robustness against complexity of the assumed CGWs emission models, as was compared in Walsh et al. (2016).
Several search methods were developed for the CGWs signals originated in isolated NSs (very extensive comparison of the sensitivity of various searches can be found in Dreissigacker et al. 2018): • The F -statistic method introduced in Jaranowski et al. (1998). The F -statistic is obtained by maximizing the likelihood function with respect to four unknown parameters of the simple CGW model of rotating NSs-CGW amplitude h 0 , initial phase of the wave Φ 0 , inclination angle of NS rotation axis with respect to the line of sight ι, and polarisation angle of the wave ψ (which are henceforth called the extrinsic parameters). This leaves a function of only four remaining parameters: f GW ,ḟ GW , δ and α (called the intrinsic parameters). Thus the dimension of the parameter space that we need to search decreases from 8 to 4. To reduce computational cost and improve method efficiency, the F -statistic can be evaluated on the 4-dimensional optimal grid of the intrinsic parameters (Pisarski & Jaranowski, 2015). As was shown in Equation (20), strength of the signal depends on the observational time: on the one hand by increasing T one can expect a detection of weaker signals, on the other hand however, analysing long-duration data requires substantial computational resources, e.g., for Polgraw time-domain F-statistic pipeline 7 , computational cost for an all-sky search scales as ∼ T 5 log(T ). Promising strategies to solve this problem are hierarchical semi-coherent methods, in which data is broken into short segments. In the first stage, each segment is analysed with the F -statistic method. In second stage, the short time segments results are combined incoherently using a certain algorithm. Several methods were proposed for the second stage: search for coincidences among candidates from short segments (Abbott et al., 2007(Abbott et al., , 2009a, stack-slide method (Brady et al., 1998;Brady & Creighton, 2000;Cutler et al., 2005), PowerFlux method (Abbott et al., 2008(Abbott et al., , 2009b with the latest significant search sensitivity improvements for O1 data (Dergachev & Papa, 2019a,b), global correlation coordinate method (Pletsch, 2008;Pletsch & Allen, 2009), Weave method Walsh et al., 2019). Independently of the details, the main goal of the F -statistic method is to find the maximum of F ( f GW ,ḟ GW , δ, α) function, and hence the parameters associated with the signal. Several optimisation procedures (such as optimal grid-based or non-derivative algorithms) were implemented in such analyses (Astone et al., 2010a;Shaltev & Prix, 2013;Sieniawska et al., 2019b). F -statistic can be evaluated on the time-domain data (Jaranowski et al., 1998;Astone et al., 2010a;Aasi et al., 2014;Pisarski & Jaranowski, 2015;Abbott et al., 2017k, 2019b and the frequency-domain data (Brady et al., 1998;Brady & Creighton, 2000;Abbott et al., 2004Abbott et al., , 2007Abbott et al., , 2009a. The main difference between these two concepts is that in the time-domain the information is distributed across the entire data set, while the frequency-domain analysis focuses on the part of the data around the frequency at which the peak appears. The data is initially calibrated in the timedomain and to be used by the frequency-domain methods, usually it is converted with the Fourier Transform methods. • The Hough transform (Hough, 1959(Hough, , 1962) is a widely used method to detect patterns in images. It can be applied to detect the CGWs signals in specific representations of the data: on the sky (Krishnan et al., 2004), and in frequencyspin-down plane (Antonucci et al., 2008;Astone et al., 2014). Both types of the Hough transform method, called SkyHough and FrequencyHough, are typically used for all-sky searches and are similar to the F -statistic are matched-filtering type methods. Due to limited computational power, they require division of data into relatively short segments. Interesting application of the Hough transform to the unknown sources searches was introduced in Miller et al. (2018). This Generalised FrequencyHough algorithm is sensitive to the braking index n, a quantity that determines the frequency behaviour of an expected signal as a function of time. In general, the evolution (decrease) of rotational frequency is described asν where K is a positive constant. Time derivative of the above equation provides the dependency of the braking index on measurable quantities (from EM observations, e.g., Espinoza et al. 2011;Hamil et al. 2015;Lasky et al. 2017;Andersson et al. 2018), frequency and its higher derivatives: Value of the braking index reveals a spin-down mechanism: n = 1 if the spindown is triggered by the relativistic particle wind; n = 3, if the spin-down is dominated by dipole radiation (as in the case of dipolar EM field); n = 5 if it is purely quadrupolar radiation (GWs emission in General Relativity); n = 7 if the spindown is due to the oscillations (lowest order r-modes, see Section 4 for further details). Some of the CGWs searches strategies assume that object is spinningdown only due to the gravitational radiation (n = 5), as it was mentioned with Equations (33) and (34), while Generalised FrequencyHough method does not assume any specific spin-down mechanism, but allows for its examination.
Hough transform has also been used in the Einstein@Home project 8 (Abbott et al., 2009a), a volunteer-based distributed computing project devoted to searching for CGWs.
• The 5-vector method (Astone et al., 2010b), in which detection of the signal is based on matching a filter to the signal + and × polarization Fourier components. The antenna response function depends on Earth sidereal angular frequency Ω ⊕ and results in a splitting of the signal power among five angular frequencies ω GW , ω GW ± Ω ⊕ and ω GW ± 2Ω ⊕ , where ω GW = 2π f GW . This method is typically used for narrowband and targeted searches.
• The Band Sampled Data (BSD) method, is dedicated for the directed searches, or those assuming limited sky regions, such as the Galactic Centre (Piccinni et al., 2019). The application of this method results in a gain in sensitivity at a fixed computational cost, as well as gain in robustness with respect to source parameter uncertainties and instrumental disturbances. From the cleaned, band-limited and down-sampled time series, collection of the overlapped short Fourier Transforms is produced. Then, the inverse Fourier Transform allows removing overlap, edge and windowing effects. Demodulation of the signal from the Doppler and spindown effects can be done e.g., by using heterodyne technique (see below). While in the F -statistic method one could manipulate with the search sensitivity by increasing the observation time, BSD method works in Fourier-domain and analogously S NR can be improved by increasing length of frequency bands (for comparison: bandwidth in F -statistic method is typically ∼0.25 Hz and in BSD ∼10 Hz).
• The time-domain heterodyne method (Dupuis & Woan, 2005) is a targeted search which uses the EM measurements of ν,ν andν (model of the phase evolution, Equation (26), assumes k = 2). The signal depends on four unknown parameters: h 0 , ψ, Φ 0 and ι. Due to the Earth's rotation, amplitude of the signal recorded by an interferometric detector is time-varying since the source moves through the antenna pattern (see Equations (8)-(10)). These variations, in the heterodyne method, are used to find characteristic frequency which is the instantaneous signal frequency, register at the detector. Additionally, frequency of the signal seen in the detector is affected due to the Earth motion. Second important step of the demodulation is to remove the Doppler shifts (correct signal time-of-arrival). A targeted search is performed with a simple Bayesian parameter estimation: first the data is heterodyned with an expected phase evolution and binned to short (e.g., 1 min) samples. Then, marginalisation over the unknown noise level is performed, assuming Gaussian and stationary noise over sufficiently long (e.g., of the order of 30 minutes) periods. 95% upper limit is defined, inferred by the analysis, in terms of a cumulative posterior, with uniform priors on orientation and strain amplitude. At the end the parameter estimation is done by numerical marginalisation. Effective and commonly used algorithm for the last marginalisation stage is called Markov Chain Monte Carlo Ashton & Prix, 2018), in which the parameter space is explored more efficiently and without spending much time in the areas with very low probability densities.
As was mentioned previously, this paper is about CGWs emission from isolated NSs. However, it is worth to mention that NSs in binary systems are also considered as a serious candidate for sources of CGWs. NSs in binary systems have an additional modulation due to the NSs movement around the binary barycenter (at least three additional parameters). To deal with the high computational cost, due to the bigger parameter space, search strategies usually rely on semi-coherent methods and are dedicated for known candidates (targeted/directed searches). Main proposed algorithms are: TwoSpect method (Goetz & Riles, 2011;Meadors et al., 2017), CrossCorr method (Whelan et al., 2015;Meadors et al., 2018), Viterbi/J-statistic (Suvorova et al., 2016(Suvorova et al., , 2017 or the Rome narrow-band method (Leaci et al., 2017). Several directed searches for CGWs from NSs in known binary systems, such as Scorpius X-1, were already performed in the past (Abbott et al., 2017o;Meadors et al., 2017).

Elastic Deformations
The simplest model of the NS CGWs emitter is described by a rigidly rotating aligned triaxial ellipsoid, radiating purely quadrupolar waves ( Figure 4). As was mentioned in Section 1.1, signal expected to be observed on the Earth will have two polarisations h × (t) and h + (t), which depend on the emission mechanism (Equation (8)): where ι is an inclination angle of NS rotation axis with respect to the line of sight, Φ(t) + Φ 0 is the phase of the wave (Φ 0 being the initial phase), expressed as a truncated Taylor series, where f (k) GW is a k-th frequency time-derivative at the Solar System Barycentre (SSB) evaluated at t = 0 ( f (1) GW =ḟ , f (2) GW =f , ...), n 0 is a constant unit vector in the direction of the NS in the SSB reference frame (it, therefore, depends on the sky position of the source) and r d is a vector joining the SSB with the detector.
The parameter h 0 is a constant GW strain, which can be estimated from Einstein's quadrupole formula (Equation (1)). For the triaxial ellipsoid model it is given by the following formula: where the deformation (also called ellipticity) is quantified by = (I 1 − I 2 )/I 3 . Quantities I 1 , I 2 , I 3 are the moments of inertia along three principal axes of the ellipsoid, with I 3 aligned with the rotation axis (see Figure 4). The symbol d denotes the distance to the source, and f GW = 2ν is the GW frequency, equal to twice the rotational frequency of the star 9 . The NS angular frequency is given by relation ω = 2πν = 2π/P, where P is the spin period. From the right-hand side of Equation (27) one can notice that the ellipticity depends on the interior properties of the NS. It carries information on how easy it is to deform NS or, equivalently, about the stiffness of the matter. So far, exact values of 9 We consider here density perturbations, which affect the spherical shape of the star δρ = Re{δρ lm (r)Y lm (θ, φ)}, where Y lm (θ, φ) denotes spherical harmonics. The multipole moment of the perturbation along radius coordinate r is Q lm = δρ lm (r)r l+2 dr. Here we focus only on the lowest-order perturbation Q 22 , consistent with l = m = 2, for which f GW = 2ν. Note that in this section we consider the simplest model, in which rotational and I 3 axes are aligned. In general they may be misaligned, producing additional CGW radiation at 1ν frequency, whose strength depends on the angle between rotational and I 3 axes and is maximal when they are perpendicular (Bonazzola & Gourgoulhon, 1996). Such cases are consider in Section 3 and 5. Searches in the LVC data for the CGW radiation at both 1ν and 2ν were performed in the past (Abbott et al., 2019g). Figure 4: Simplest NS CGW model: non-axisymmetric, rotating NS (described as a triaxial ellipsoid) radiating gravitational waves at twice the spin frequency.
are uncertain. The empirical formula for the moment of inertia from Bejger & Haensel (2002), obtained for a set of realistic EOS, is where the compactness is redefined as x = (M/M )(km/R), and a(x) functions are characteristic for the hadronic NSs and strange stars (SSs, composed of deconfined quarks, discussed below in this section): According to Owen (2005), maximum ellipticity NS max depends on breaking strain of the crust σ max , mass M and radius R in the following manner: (30) For the typical NS parameters (M = 1.4 M , R = 10 km), the maximum ellipticity can be approximated as : However, σ max values are also ambiguous: if the NS crust is a perfect crystal with no defects, the σ max ≈ 0.1 (Kittel, 2005;Horowitz & Kadau, 2009;Chugunov & Horowitz, 2010). For more amorphous materials, the breaking strain is expected to be much smaller: σ max ≈ 10 −5 − 10 −2 (Ruderman, 1992;Kittel, 2005). Nevertheless, Equation (31) suggests that even for the extreme value of σ max ≈ 0.1, one can expect NS max 5 × 10 −6 .
So far, we have assumed that NSs are made of hadronic matter, whose EOS is largely unknown. Another hypothesis about the dense NS matter was proposed in the literature almost fifty years ago and is still waiting to be proven: Ivanenko & Kurdgelaidze (1965); Itoh (1970);Fritzsch et al. (1973); Baym & Chin (1976); Keister & Kisslinger (1976); Chapline & Nauenberg (1977); Fechner & Joss (1978) considered a possibility that in extreme conditions, e.g., under enormous pressure P, baryon matter is unstable. Protons and neutrons are made of two types of quarks: u-up and d -down (proton = uud; neutron=udd) and forces between quarks are mediated by gluons. It is intuitive that the mixture of protons and neutrons can deconfine only into 2-flavour quark matter (baryons deconfine into quasi-free u and d quarks). Hypothetical objects made of such matter are now called quark stars (QS).
However, when chemical potential grows, about half of d quarks can be transformed into s-strange quarks. Such conversion is due to the weak-interaction process, in contrast to the conversion of baryons into quarks, which is a strong-interaction process (see Chapter 8 in Haensel et al. 2007 for the review). Special attention will be paid to the sub-class of QS, called strange stars (SSs), in which strange quarks-s coexist with quarks d and u (Bodmer, 1971;Witten, 1984;Haensel et al., 1986;Madsen, 1998). This unique 3-flavour mixture ( which is in weak-interaction equilibrium) would be more stable than 2-flavour quark matter, but also more stable than hadronic matter such as, e.g., iron 56 Fe (which is one of the most tightly bound nuclei). For the detailed information see a review by Weber (2005). Some authors hypothesize that SSs may be in a solid state (see e.g., Xu 2003), having also a solid crust made of 'normal' matter. Contribution of such crust to the total ellipticity will be order of a few percent correction to that of the internal layers.
Similarly to Equation (30) From the above, one sees that S S max is expected to be a few orders of magnitudes larger than NS max . Broad study on maximal ellipticity for multiple EOS for the Newtonian and relativistic stars was performed in Haskell et al. (2007); Mannarelli et al. (2007Mannarelli et al. ( , 2008; Knippel & Sedrakian (2009);Johnson-McDaniel & Owen (2013), including superconducting quark matter (SQM) for which S QM max can be smaller than 10 −1 . Observational determination of max value would put strong constraints on the type and properties of the dense matter inside compact objects.
For more than 200 of known pulsars with known frequency and frequency derivative from EM observations, one can make the following assumption to illustrate the order of magnitude. Let us assume that the observed pulsars are spinning down solely due to loss of the gravitational energy. Such a hypothetical object is denoted as the gravitar (Knispel & Allen, 2008). For a given detector sensitivity, an indirect spindown limit can be thus established (Jaranowski et al., 1998): Combination of Equations (27) and (33) results in a corresponding spin-down limit for the fiducial equatorial ellipticity: So far, in the O1 and O2 LVC observing runs, no CGWs signals were detected, but instead interesting upper limits in a broad range of frequencies for known pulsars were set (Abbott et al., 2017n, 2019c. For example, the currently best limits for Crab pulsar (J0534+2200) are h sd = (1.4 ± 0.4) × 10 −24 resulting in sd = 7.56 × 10 −4 . Additionally, in the O2 data analysis, for some of the pulsars (J1400-6325, J1813-1246, J1833-1034, J1952+3252, J0940-5428, J1747-2809, J1400-6325, J1833-1034, J1747-2809), the energy spin-down limits (the maximal available rotational energy loss due to CGW emission) were surpassed (Abbott et al., 2019c). For example, the Crab pulsar (J0534+2200) is emitting ∼1% or less energy in CGW (for more details see Table IV in Abbott et al. (2019c)).
Prospective detection of CGWs, consistent with a deformed NSs, will be a breakthrough in several topics, mostly related to the properties of matter under huge pressures. It will put strong constraints on EOS, and likely may yield a distinction between compact stars built of the hadronic and strange matter. It will also allow for testing of the outer NS layer-the crust. An analysis of the crustal failure and its other properties may reveal not only astrophysically important information, but also can be valuable in engineering, e.g., mechanics and material strength studies. Finally, it allows measuring deformation itself, describe the exact NS shape and mechanisms that triggered deviations from a spherical shape.

Magnetic Field
In a widely accepted pulsar model, inferred from the EM observations, misalignment between the global dipole magnetic field axis and the rotation axis is responsible for observed pulsations. It is likely that for some range of conditions (described in this section), an asymmetry in the magnetic field distribution in the interior of NS may lead to the emission of CGWs.
The idea that magnetic stresses can deform a star and lead to the CGWs emission was originally proposed by Chandrasekhar & Fermi (1953) (1984), where a NS was approximated by a rigidly rotating Newtonian incompressible fluid body, with the uniform internal magnetic field and the dipole external magnetic field. Inclination of the magnetic dipole moment with respect to the rotation axis is given in this model by angle χ. Note that for misaligned rotational and magnetic axes, time-varying quadrupole moment will also be present for I 1 = I 2 . In such conditions, the NS is a biaxial ellipsoid, illustrated in Figure 5. The ellipticity is then described by the ratio = (I 3 − I 1 )/I 1 , and the resulting CGW strain will take form (Bonazzola & Gourgoulhon, 1996): where coefficient β measures efficiency of the magnetic structure in distorting the star (magnetic distortion factor). For a simplified model (incompressible fluid, uniform internal magnetic field), parameter β equals 1/5. Information about the induced deformation is encoded in β factor, as = βM 2 /M 2 0 , where M is the magnetic dipole moment and M 0 has the dimension of a magnetic dipole moment in order to make β dimensionless. It is clear that h 0 grows when χ → 0 (but for practical reasons χ cannot be too small because it will break the formula down), and/or for large β. As was shown in Bonazzola & Gourgoulhon (1996), such emission occurs mostly on f GW = 1ν. One can also notice, by comparison of the prefactors of Equations (27) and (35) that the expected CGW strain from magnetic effects is almost four orders of magnitude smaller than in the case of triaxial ellipsoid. For the above model, in case of the incompressible NS with poloidal magnetic field ellipticity is expressed as (Haskell et al., 2008): whereB is the volume averaged magnetic field and factor A is a constant, characteristic to the model (e.g., A = 10 −12 in case of a uniform density star, A = 8 × 10 −11 for the polytrope n = 1 with the purely poloidal field, A = −5 × 10 −12 for the polytrope n = 1 with purely toroidal field, etc.) 10 . It was shown that purely poloidal or purely toroidal magnetic fields are unstable (Wright, 1973;Braithwaite & Spruit, 2006;Braithwaite, 2007), and that the realistic description of NS requires mixed field configurations (Braithwaite & Nordlund, 2006;Ciolfi et al., 2009;Lander & Jones, 2012), giving an estimate of ellipticity (Mastrano et al., 2011): where B poloidal is the poloidal component of the surface magnetic field, and Ξ is the ratio of poloidal-to-total magnetic field energy: Ξ = 0 corresponds to purely toroidal field and Ξ = 1 to purely poloidal field. NSs have magnetic fields below 10 15 G (Shapiro & Teukolsky, 1983), so the effect of magnetic deformation is significantly smaller than the one discussed in Section 2. NSs with magnetic fields typically of the order of 10 15 G are called magnetars. Known magnetars rotate with periods of P∼ 1-10 s (Manchester et al., 2005), so the frequency of the emitted CGWs would be very low, in the frequency range unavailable for Advanced LIGO, Advanced Virgo and next generation detectors, such as the Einstein Telescope (Lasky, 2015). Ref. Cutler (2002) and recent evaluations for different EOS and bulk viscosity models in Dall'Osso et al. (2018) show that, while the first-year spin-down of a newborn Figure 5: A model of the CGWs emission due to the magnetic deformation. The angle χ measures the misalignment between the magnetic field axis and rotation axis. In principle for χ 0, CGWs emission will be produced even if I 1 = I 2 . Note that on this schematic picture B represents the dominant toroidal field, reflecting the fact that NS shape is prolate along the magnetic axis.
NS is most likely dominated by EM processes, reasonable values of internal toroidal field B toroidal and the external dipolar field B external can lead to detectable GWs emission, for an object in our Galaxy. While the centrifugal force distorts a NS into anoblate shape, the internal toroidal magnetic field makes them prolate. B toroidal determines the NS shape if B toroidal 3.4 × 10 12 ν 300 Hz 2 G. (38) The above formula is important not only for the newborn NSs, but also for the NSs in binary systems, where accretion processes are active and affecting magnetic field. That may happened e.g., in low-mass X-ray binaries (LMXB)-systems in which mass transfer between the companion and NS is via Roche-lobe overflow. Material in accretion disc formed around NS is heated so much that it shines brightly in X-rays. NSs in LMXB initially possess B external ∼ 10 12 -10 14 G, B toroidal ∼ 10 12 -10 15 G and their spin and magnetic axes are nearly aligned (θ ≈ 0), as was studied in Cutler (2002). Such objects spin-down electromagnetically, until accretion process from the companion is ongoing: accretion reduces B external below 10 9 G. However, in the interior, B toroidal remains unchanged. Dissipation processes dominate and at some point magnetic axis rapidly 'flips' orthogonally to the rotational axis. For the full understanding of the LMXB evolution and GWs emission, one should consider comparison between timescales of four processes: (i) spin-down due to the GWs emission: whereṀ is the accretion rate; (iiii) dissipation timescale on which the instability acts: where Y = τ DIS S B P is a parameter, which was found hard to estimate and B is the deformation caused by the toroidal field and defined as: As it was shown in Dall'Osso et al. (2018), for an optimal range of B ∼ (1−5)×10 −3 and birth spin period 2 ms the horizon of CGWs detectability is equal to 4, and 30 Mpc, for Advanced and third generation interferometers at design sensitivity, respectively. Additionally, LMXB give opportunity to estimate maximum expected signal from accreting neutron stars, taking advantage from EM observations. If the angular momentum lost in GWs is recovered by accretion, then the strongest GWs emitters are those accreting at the highest rate, namely LMXB. As it was shown in Papaloizou & Pringle (1978); Wagoner (1984); Bildsten (1998), if the spin-down torque from GWs emission is in equilibrium with the accretion torque, h 0 is directly related to the X-ray luminosity (Wagoner, 1984;Bildsten, 1998;Ushomirsky et al., 2000): where F x is the observed X-ray flux; note that here the information about the distance is encoded in F x and both-GWs and X-rays-are falling inversely proportional to d 2 . Most of LMXB have spins in the relatively narrow range 270 Hz ν 620 Hz (Chakrabarty et al., 2003), while the Keplerian frequency (discussed in Section 1.4) is typically much bigger, ∼1400 Hz. That leads to the conclusion that LMXB have to be equipped in an additional mechanism balancing the angular momentum transfer and preventing spin increase. According to Equation (44), brightest known X-rays source, Sco X-1 (Giacconi et al., 1962) should also emit the strongest GW signal: Unfortunately, Sco X-1 spin frequency is not well constrained, it is nevertheless one of the prime targets for LVC and its GW signal was searched for in the past (Abbott et al., 2007;Aasi et al., 2015a;Abbott et al., 2019f). During the latest and the most sensitive search (Abbott et al., 2019f), no evidence for the GWs emission was found; 95% confidence upper limits were set at h 95% 0 = 3.47 × 10 −25 , assuming the marginalisation over the source inclination angle.
As was mentioned several times in this section, CGWs created by the magnetic processes typically have amplitudes smaller than the sensitivity of Advanced LIGO and Advanced Virgo. However, it may be within the reach of a future improved network of detectors. Successful detection of such signals would be an amazing tool for probing the magnetic fields in NSs: their composition (poloidal, toroidal), strength and evolution. Indirectly one could test deformability and compressibility of the NSs. These studies will be also complementary to EM observations of special NSs classes: magnetars, young NSs and NSs in LMXB systems. Closer look through the GW analysis could possibly allow for testing of accretion, precession and spin evolution processes.

Oscillations
In asteroseismology, the Rossby modes (also called r-modes, Rossby 1939) are a subset of inertial waves caused by the Coriolis force acting as restoring force along the surface. In the NSs the r-modes are triggered by the Chandrasekhar-Friedman-Schutz instability (Chandrasekhar, 1970;Friedman & Schutz, 1975). This instability is driven by GW back-reaction-it creates and maintains hydrodynamic waves in the fluid components, which propagate in the opposite direction to that of the NS rotation, producing GWs. So far, the EM observations of the oscillations in NSs are insufficient in the context of EOS, because it is impossible to directly observe the modes or thermal radiation from the surface. Details of the methods and limitations of EM asteroseismological observations can be found e.g., in Cunha et al. (2007). It is expected that the CGWs observations will complement our understanding of the compact objects oscillations.
The r-modes were first proposed as a source of GWs from newborn NSs  and from accreting NSs (Bildsten, 1998;Andersson et al., 1999). It was shown that the r-modes can survive only in a specific temperature window, in which they remain unstable: at too low temperatures, dissipation due to shear viscosity damps the mode and when the matter is too hot, bulk viscosity will prevent the mode from growing, as shown in the schematic plot in Figure 6. R-mode instability window is open for NSs with temperatures between 10 6 -10 10 K Bildsten, 1998;Bondarescu et al., 2007;Haskell, 2015). This information allows testing properties of the interiors of NSs and theoretical models that describe matter interactions in such conditions. The simplest model 11 of the CGW emission, triggered by the unstable rmodes in the newborn NSs, was introduced in Owen et al. (1998), where authors shown that CGW amplitude in this case can be expressed by: where α is dimensionless r-mode amplitude andρ is the mean density.
Also in the temperature range that corresponds to the r-modes instability window, conditions to form superfluidity are fulfilled. The superfluid component plays an important role in NSs interiors e.g., in attempts to explain glitches e.g., sudden increases of rotational frequency (these events are regularly observed, but not explained theoretically, Haskell 2015). For the NSs in LMXB systems, the r-modes can affect the spin evolution of rapidly spinning objects (Haskell, 2015;. This mechanism was proposed not only for young and hot NSs, but also for the old, accreting NSs. In the LMXB, together with matter also the angular momentum is 'accreted', resulting in spinning up the object. This is thought to be the mechanism by which old, long-period NSs are recycled to millisecond periods (Alpar et al., 1982). In principle, long-lasting persistent accretion may spin NSs up to the Keplerian frequency (see previous Section), but a spin frequency cut-off of about 700 Hz is observed (Chakrabarty et al., 2003). One of the proposed mechanisms to balance the spin evolution is the GWs emission due to the r-mode instability (Bildsten, 1998;Andersson et al., 1999). When the r-mode amplitude, defined in Equation (53), is very large, the system will be located deeply inside the instability window, but also the evolution will be very rapid and the system will rapidly evolve out of the window (Levin , 1999;Spruit, 1999). On the contrary, when the r-mode amplitude is low, the evolution is slower, but system will remain close to the instability curve (Heyl, 2002). As a consequence, in both scenarios it is highly unlikely to observe a LMXB deeply inside the instability window. However, according to Haskell (2015), our theoretical understanding is not fully consistent with observations: a large number of observed LMXBs is located deeply in the r-mode instability window (Haskell et al., 2012). Mechanisms, which can potentially explain this inconsistency, are e.g., the presence of exotic particles in the core (Alford et al., 2012;Haskell et al., 2012), viscous damping at the crust/core boundary layer interface (Glampedakis & Andersson, 2006;Ho et al., 2011) or strong superfluid mutual friction Haskell et al., 2009). An advantage of the LMXB systems is that they are observed in the EM spectrum and some systems were tentatively interpreted as sources of r-modes recently, e.g., by Andersson et al. (2018). The study of r-modes provides a huge opportunity to test various models of NSs interiors and the physics behind their spin evolution.
For the CGWs triggered by r-modes oscillations, the (Newtonian) relation for the lowest order (strongest) mode 12 is f GW = 4/3ν. This relation is independent on EOS and it is a good approximation for the slowly rotating NSs. However, according to numerical results, inclusion of relativistic corrections may increase f GW by a few tens of percent (Lockitch et al., 2003;Pons et al., 2005;Idrisy et al., 2015;Jasiulek & Chirenti, 2017). Exact value depends mostly on the object compactness C (dimensionless massradius ratio: C ≈ 0.207 M 1.4 M 10 km R , for NSs typically 0.11 C 0.31), but also on rotational frequency of the object (Lindblom et al., 1999;Idrisy et al., 2015;Jasiulek & Chirenti, 2017), including rigid versus differential rotation (Stavridis et al., 2007;Idrisy 12 Unlike in the case of the elastic deformations, what was shown in Section 2, r-modes affect equation of motion δv = α(r/R) l RωY B lm exp(iω r t), where α is dimensionless amplitude and ω r = − 2m l(l+1) ω is the angular frequency in co-rotating frame. The only non-trivial solution is for l = m = 2, giving ω r = − 2 3 ω, what can be transferred to the angular frequency ω i in the inertial frame: ω i = − 2 3 ω + 2ω = 4 3 ω. Dividing this expression by 2π gives f GW = 4 3 ν.  Levin & Ushomirsky, 2001) and EOS (Lockitch et al., 2003;Pons et al., 2005;Lattimer & Prakash, 2007;Idrisy et al., 2015). Additionally, as it was shown in Arras et al. (2003), consideration of the non-linear coupling forces among internal modes leads to the result that r-mode signal from both newly born NSs and LMXB in the spin-down phase of Levin's limit cycle 13 (Levin , 1999) will be detectable by enhanced LIGO detectors out to ∼100-200 kpc. Dissimilarities of the non-linear effects in supra-thermal regime for hadronic and quark EOS were studied by Alford et al. (2010). Nonetheless, for practical purposes (such as CGWs searches), a Newtonian approximation is typically used. Potentially, the detection of r-modes can confirm or rule out the existence of SSs, introduced and discussed in Section 2. The viscosity coefficients of NSs and SSs differ significantly and as a consequence, the r-modes in a SS are unstable at lower temperatures (Chatterjee & Bandyopadhyay, 2008;Alford et al., 2010;Haskell et al., 2012). The existence, evolution and properties of SSs is still an open question, widely discussed in the literature (see e.g., Jaikumar et al. 2006;Blaschke et al. 2002). As shown in Mytidis et al. (2015), a hypothetical r-mode detection can put constraints on the moment of inertia of the compact object and bring closer to understanding the EOS.
The growth time of the instability for r-modes can be relatively short 14 (Andersson & Kokkotas, 2001;Haskell, 2015): while the characteristic damping timescales associated with the bulk viscosity, τ bv , and shear viscosity, τ sv are where T is the NS core temperature. In this minimal model it is assumed that at high temperatures, bulk viscosity due to modified URCA reactions provides the main damping mechanism, while at low temperatures the main contribution is from shear viscosity, due to electron-electron scattering processes. Instability curve from Figure 6 can be calculated with the following simple formula: where 1/τ diss = 1/τ bv + 1/τ sv + additional processes. Several additional mechanisms were considered in the literature, for example the crust/core velocity difference (Levin & Ushomirsky, 2001;Glampedakis & Andersson, 2006), exotic particles in the core (Andersson et al., 2010;Alford et al., 2012), or strong superfluid mutual friction Haskell et al., 2014). These models can lead to very significant changes in the shape, width and depth of the instability window. Another natural and astrophysically motivated targets are the binary systems in which the accretion is present. CGW strain from the r-modes excited in accreting systems can be estimated as (Owen, 2010;Chugunov, 2019): GWs emitted by a NS destabilised by r-modes carry the information about the internal structure and physical phenomena inside the star and are, in principle, a very powerful tool for testing extreme conditions of the NSs interiors, providing information that cannot be obtained by other means. Constraints on the EOS were considered plest illustrative model. General expressions of the timescales are is the shear viscosity factor; 1 , where ζ = 6.0 × 10 −59 3 2ω 2 ρ 2 T 6 is the bulk viscosity factor. In principle all these timescales are sensitive to EOS due to the occurrence of density ρ in the equations, by e.g., Mytidis et al. (2015); the superfluid layer and crust-core interface properties by Haskell (2015); rotational frequency evolution in accreting binaries by Andersson et al. (1999), spin-down of the young NSs by Alford & Schwenzer (2014). Similarly to Section 2, one can define the spin-down limit for r-modes, obtained by assuming that all of the observed change in spin frequency is due to the GWs emission (Owen, 2010): 1.6 × 10 −24 1 kpc d |ḟ GW | 10 −10 Hz/s 1/2 100 Hz and the corresponding r-mode amplitude parameter : From the point of view of detecting the r-modes, several methods were used so far to look for their signatures in the LVC data. For example, the F -statistic method described in Section 1.5 is so general that can be applied also for r-modes CGWs searches, e.g., during supernovae remnant searches (Abbott et al., 2019h). No signal was found, but interesting fiducial r-mode amplitudes upper limits (α sd ≈ 3 × 10 −8 ) were set.
Searches for r-modes from pulsars with known sky position require at least three free parameters f GW ,ḟ GW ,f GW . Fortunately, for some pulsars these parameters are well measured. As shown in Caride et al. (2019), the selection of appropriate ranges of frequencies and spin-down parameters is crucial. In such a case, for most pulsars from ATNF Pulsar Database, number of required GWs templates is ∼ 10 9 -10 12 , comparable in terms of the computational cost to the CGWs searches performed in the past in the O1 and O2 observational runs for triaxial ellipsoid models (see Section 2).
Additionally, a special type of oscillations may come from the newborn NSs. Following the gravitational collapse, the proto-neutron star (PNS) radiates its binding energy (about 0.1 M ) via neutrino emission in a timescale of tens of seconds. A small fraction of this energy can be released through violent oscillations leading to GWs emission. As shown in Ferrari et al. (2003), initial amplitude of such oscillations for the Galactic PNSs should be larger than h 0 10 −22 that it is detectable by the LVC, what is marginally consistent with numerical studies of the axisymmetric collapse of the core of a massive star (Dimmelmeier et al., 2002). However, such a signal is expected to last for a very brief moment (tens of seconds) and its waveform will depend on multiple variables, such as the poorly known high-temperature dense-matter EOS, resulting with a computationally expensive and uncertain data analysis results. Additionally, Galactic supernovae events are rare (approximately one per century), making the PNSs GWs unlikely in the LVC data.
It is worth metioning that r-modes are not the only oscillation type that may, in principle, lead to the CGWs emission. Initially, the Chandrasekhar-Friedman-Schutz instability was considered a source of the fundamental fluid mode: the f-mode. However, to produce CGWs from the f-mode instability, very fast rotational frequency is needed, around 95% of the Keplerian frequency (Friedman, 1983;Ipser & Lindblom, 1991), as well as hot, T 5 × 10 8 K in accordance with recent models for the observed cooling of Cassiopeia A (Page et al., 2011;Shternin et al., 2011), non-superfluid matter (Ipser & Lindblom, 1991;Passamonti et al., 2013). Such a situation could be possible only if superfluid core is not formed yet, which would correspond to rapidly spinning, newborn NSs. Nevertheless, further works by Gaertig et al. (2011);Doneva et al. (2013) unveil that relativistic, massive (M ≈ 2 M ) NSs with realistic EOS can support a wider instability window than it was initially thought, especially for the l = m = 4 multipole. Another consequence of that result is crucial in the context of binary NS mergers: components of binary NS coalescence and also final product of the merger heat up. If the remnant is supermassive and fast-spinning (close to its Keplerian frequency) unstable NSs, it may emit post-merger CGW signal. For the reviews on NSs oscillations types, see Kokkotas (1997); Glampedakis & Gualtieri (2018).
CGWs from oscillations will carry asteroseismological information about NSs interior: its hydrodynamics, superfluidity properties, viscosity and temperature. They may indirectly allow for a distinction between NSs and SSs. Detection of such signals will be invaluable in understanding of the details of newborn and accreting NSs.

Free Precession
Historically, freely precessing NSs were first considered as good sources of CGWs emission by Zimmermann (1978); Alpar & Pines (1985). Let us adopt a free precessing NS model consisting of a shell containing a liquid. For a perfectly elastic shell and the inviscid fluid, there will be no energy dissipation: such star will precess with a constant angular velocity and wobble angle.
However, a realistic model requires realistic elasticity and fluidity description, as well as the internal energy dissipation and external, astrophysical torque mechanisms balancing energy dissipation. Internal dissipation means here that energy losses inside NS are converting the mechanical energy into heat. Strong interactions between vortices and the normal component (mutual friction, mentioned already in Section 4), are efficient dissipative channels that can dumped precession (Shaham, 1977;Sedrakian et al., 1999;Haskell & Sedrakian, 2018b). Reasoning for a mutual friction in the case of type-I superconducting protons is difficult because of lack of a model-independent predictions for the domain structure and size of type-I superconductor (Haskell & Sedrakian, 2018b). However, as it was shown in Wasserman (2003), if the magnetic stresses are large enough, precession is inevitable. Such enormous magnetic stresses can arise if the core is a type-II superconductor or from toroidal fields ∼ 10 14 G, if the core is a normal conductor (Wasserman, 2003). Observational evidence of free precession can put strong constraints for the mutual friction and confirm weak coupling of the superfluid to the normal component.
Additionally, if the magnetic and rotational axes are misaligned, NS will precess (Mestel & Takhar, 1972). That process can be especially important in the early life of the millisecond magnetars (Lander & Jones, 2017. For the very young NSs, before the crust solidifies, any elastic component in the moment of inertia tensor is allowed and effective ellipticity 15 B (introduced in Section 3, Equation (43)) comes entirely from the magnetic deformation. When considering precession, χ from Figure 5, between axis of rotation and I 3 becomes the angle around which the rotational axis moves around. It is commonly called the wobble angle and denoted by θ in Figure 7. Rotational frequency of the slow precession Ω p (assuming rigid rotation of the uniform-density object) can be expressed as (Mestel & Takhar, 1972): where ω = 2πν. For non-rigid rotation, the centrifugal force can deform a star by itself; a quantity corresponding to size of a centrifugal distortion is α and both deformation components are related to the global parameters of the object (Lander & Jones, 2017: Dynamical evolution of the magnetic field in the newborn NSs is crucial to understand their GWs emission. Interplay between energy loss due to the inclination angle evolutionθ and the internal energy-dissipation rate due to viscous processes (see Section 3) may lead to the precession damping and change of angle θ. As it can be deduced from Equation (35), maximal CGWs emission can be expected when rotational and magnetic axes are orthogonal to each other. Ref. Lander & Jones (2018) shows that below 10 14 G, all NSs at some point of evolution have orthogonal rotational and magnetic field axes, regardless of their birth spin frequency. Above 10 14 G only those NSs that are born spinning fast enough can enter the orthogonalisation region. Additionally, energy and angular momentum are carried away as GWs from a freely precessing NSs (Bertotti & Anile, 1973). If the body is rigid enough, θ diminishes monotonically, but elastic behaviour may instead increase θ to π/2. While dissipation mechanisms tend to decrease θ, pumping mechanisms, such as e.g., accretion, are increasing θ. As shown in Lamb et al. (1975) accretion torque can be effective in exciting large amplitude precession of a rigid body with the excitation timescale (τ ACC , introduced in Section 3) comparable to the spin-up timescale. The criterion for excitation is τ excitation < τ damping . The evolution of the wobble angle can be therefore expressed asθ where in the simplest model one can assume τ excitation = τ ACC and τ damping = τ CP . In principle, inequality τ excitation < τ damping can be used to set upper limits for the GWs amplitudes (for the NSs with known or assumed rotational frequency, accretion rate and crustal breaking strain). Unfortunately, such a signal will be a few orders of magnitude too weak to be detected by the LVC and the Einstein Telescope (Jones & Andersson, 2002).
inertia of the spherically symmetric density field ρ 0 . GW strain for free precession is given by Zimmermann (1978); Jones & Andersson (2002); Van Den Broeck (2005): where ∆I d = I 3 − I 1 is the strain-induced distortion of the whole star, also called effective oblateness. Wobble angle θ cannot be too big, because too big a deformation of the object (too large breaking strain σ max introduced in Section 2) may destroy NS, so θ max ∼ σ max . Generally, the CGWs emission originating in free precession process is present at frequencies f GW,1 = 2ν, f GW,2 = ν + ν p , where ν p is the precession frequency (Zimmermann, 1978); after including second order expansion (∼ θ 2 ) also f GW,3 = 2(ν + ν p ), according to Van Den Broeck (2005). This third frequency will be seen in GWs spectra with h 0 almost two orders of magnitude smaller than first and second frequencies. According to Van Den Broeck (2005), f GW,1 should be detectable if NS spin-down age is much less than 10 3 yr and observability of f GW,2 and f GW,3 depends on the crustal breaking strain, which is very uncertain parameter, see Section 2.
To extract information about the wobble angle θ and the deformation at least two of the above-mentioned, characteristic frequencies need to be detected. Free precession may be much longer lived (∼ 10 5 years) than initially thought (Cutler & Jones, 2000), resulting with weak CGWs emission h 0 ∼ 10 −28 or smaller. According to Jones & Andersson (2002), it is impossible to find astrophysical pumping mechanisms capable of producing CGWs detectable by an Advanced LIGO.
Some hints about the precession were delivered by EM observations of the radio pulsars. The most convincing evidence for free precession is provided by 13 years observations of PSR B1828-11. Variations of arrival-time residuals from PSR B1828-11 were interpreted as a precession, with precession period ≈ 250 days (Stairs et al., 2000;Akgün et al., 2006). Results obtained from observational data are consistent with the model of precession of a triaxial rigid body, with a slight statistical preference for a prolate NS shape (Akgün et al., 2006). Also analysis of the timing data from another source, PSR B1642-03, also were interpreted as a precessing NS, with characteristic variation from 3 to 7 years (Cordes , 1993;Shabanova et al., 2001).
Free precession may cause deformations of the object, subsequently allowing testing the NS dense-matter properties, such as breaking strain, viscosity, rigidity and elasticity. CGW detection from free precession will allow also for better understanding of processes present inside the NSs, such as heating, superfluidity, dissipation, mechanisms of the torque and its evolution.

Summary
In recent years we have witnessed the first direct detections of gravitational waves (Abbott et al., 2016b(Abbott et al., , 2017a(Abbott et al., ,b,c,d, 2018b. It started a new observational field: the gravitational wave astronomy. Those detections corresponded to observation of coalescences of the binary systems: violent and energetic phenomena involving black-holes, impossible to detect with traditional astronomical methods, as well as electromagnetically bright binary neutron stars system coalescence. As the interferometers improve their sensitivity (Harry et al., 2010;Acernese et al., 2014;Moore et al., 2015), we expect other types of signals to be detected. One promising scenario is the detection of continuously emitted, periodic and almost-monochromatic GWs produced by rotating NSs. Several mechanisms can be responsible for such a GWs emission: elastically and/or magnetically driven deformations (mountains on the NS surface supported by the elastic strain or magnetic field), thermal asymmetries due to accretion, free precession or instabilities leading to modes of oscillation (r-modes).
Discovery of a persistent source will be the capstone of GW astronomy will allow for repeatable observations, which will result in acquiring more knowledge and better understanding of the NSs interiors, especially on their crust physics and elastic properties at low temperatures (study of cold EOS as opposed to hot EOS in binary NS collisions). In general, these NSs may have different masses and matter conditions than the objects in binary NS inspirals and mergers, which leads to valuable GWs observations in different dynamical regimes and, additionally, to various tests of the general relativity by exploiting the fact that detectors are moving with respect to the source (e.g., studies of independent wave polarizations). In principle, NSs observed in CGWs could be also a valuable tool in the detectors' strain calibration, astrophysical distance ladder calibration (useful in cosmology) and measurements of fundamental aspects of gravity Isi et al., 2017).
The most promising CGWs emission scenario (giving the highest GW amplitude) assumes deformations on the NS surface. Because of that, most of the LVC effort is concentrated on searches of signals consistent with the model of triaxial ellipsoid radiating CGWs at twice the spin frequency (Abbott et al., 2017k, 2018c(Abbott et al., 2017k, , 2019b(Abbott et al., 2017k, , 2017c. However, in the future, when LIGO and Virgo detectors will reach their design sensitivity and new instruments (ground-and space-based, Akatsu et al. 2017;Amaro-Seoane et al. 2017;Sathyaprakash et al. 2012;Unnikrishnan 2013) will join the network, testing of other models and CGWs sources will be possible.
Additionally, by combining CGWs and EM observations one can break the degeneracy in some expressions considered in this paper. For example, for the triaxial ellipsoid model, by knowing distance and spin frequency from EM observations and CGW amplitude from GWs detections, immediately one can obtain size of the deformation. This quantity is a few order of magnitudes different for NSs and SSs and carries information about the crust properties. For the strongly magnetised objects, such as magnetars, newborn NSs or NSs in LMXB systems joint EM and CGWs analysis will allow for testing magnetic field composition, strength and evolution. Also unstable r-modes carry information about the NSs interiors: their superfluid layer, viscosity and temperature. Finally, the precession-in the case of detection may expose information about magnetic field, elasticity and superfluidity of the NSs.
Potential CGWs detection will allow for testing NSs matter (their EOS, crustal properties, deformability), processes inside the object (heating, superfluidity, instabilities, magnetic field evolution), close environment of the stars (their magnetic field, accretion processes). Such unique analysis will be essential also in our understanding of NSs evolution, especially for the newborn objects. Very young NSs undergo several dynamical mechanisms and usually are unstable -such changes are difficult to observe in traditional, EM telescopes. GW astronomy will open a new door for NSs studies and may be a crucial piece in solving the mystery of NSs nature and EOS.

Funding
The work was partially supported by the Polish National Science Centre grants no.