The Gravitational-Wave Physics

The direct detection of gravitational wave by Laser Interferometer Gravitational-Wave Observatory indicates the coming of the era of gravitational-wave astronomy and gravitational-wave cosmology. It is expected that more and more gravitational-wave events will be detected by currently existing and planned gravitational-wave detectors. The gravitational waves open a new window to explore the Universe and various mysteries will be disclosed through the gravitational-wave detection, combined with other cosmological probes. The gravitational-wave physics is not only related to gravitation theory, but also is closely tied to fundamental physics, cosmology and astrophysics. In this review article, three kinds of sources of gravitational waves and relevant physics will be discussed, namely gravitational waves produced during the inflation and preheating phases of the Universe, the gravitational waves produced during the first-order phase transition as the Universe cools down and the gravitational waves from the three phases: inspiral, merger and ringdown of a compact binary system, respectively. We will also discuss the gravitational waves as a standard siren to explore the evolution of the Universe.

produced in the primordial Universe and its detection through CMB. In particular, we emphasize the properties of GWs created during inflation and preheating processes and current situation of detection of the primordial GWs. In section 3 we will discuss the GWs produced during cosmic PT. There we will introduce the bubble nucleation, bubble expansion and bubble percolation. During the strong first-order PT, bubble collision, turbulent magnetohydrodynamics (MHD) and sound wave are all the sources to produce GWs. Section 4 will mainly be devoted to discussions on the GWs from the dynamics of compact binary systems. There we will introduce three main methods to solve the binary system: the post-Newtonian (PN) approximation, numerical relativity and the black hole perturbation, corresponding three phases: inspiral, merger and ringdown of two black holes system, respectively. In that section, we will also discuss the possibility of GWs produced by binary systems as a cosmological probe. With this new probe, it is expected to have a strong constraint on cosmological parameters, combining with other cosmological probes.

Gravitational Waves From Primordial Universe
In order to solve some problems in big-bang cosmology such as the horizon and flatness problems, inflationary scenario was introduced [8][9][10], in which a period of accelerated expansion of the Universe happened at early times. Inflation not only predicts the primordial scalar perturbations, which provide a natural way to generating the anisotropies of the CMB radiation and the initial tiny seeds of the large-scale structure observed today in the Universe, but also generates a stochastic background of the primordial GWs. Although such a stochastic background of the primordial GWs has not been observed yet, its detection would open a new window to understanding the physics of the early Universe and thus the origin and evolution of the Universe. In this section, we shall firstly review the properties of the primordial GWs produced during inflation and preheating (see Fig. 1), and then discuss observational implications.
GWs are described by a transverse-traceless gauge-invariant tensor perturbation, h ij , in a Friedman-Robertson-Walker (FRW) metric, ds 2 = a 2 (τ ) −dτ 2 + (δ ij + h ij )dx i dx j , (2.1) where τ is the conformal time and a is the scale factor, and h ij satisfies ∂ i h ij = 0 and δ ij h ij = 0. To first order in h ij , the perturbed Einstein equation reads where the prime denotes the derivative with respect to τ , H ≡ a /a is the Hubble parameter in τ , M pl ≡ (8πG) −1/2 is the reduced Planck mass, and the source term; Π TT ij is the transverse-traceless projection of the anisotropic stress tensor T ij . Since h ij is symmetric, transverse and trace-free, tensor modes are left with two physical degrees of freedom, which are expanded in the Fourier space as ϕ Figure 1. The schematic illustration of the inflaton potential. In the slow-roll inflationary scenario, an accelerated expansion of the Universe occurs when the inflaton rolls slowly along its potential. After inflation ends, the inflaton oscillates around the minimum of its potential, whose energy is converted into radiations.

During Inflation
In the standard single-field slow-roll inflationary scenario, at first order in perturbation theory Eq. (2.2) reduces to a free wave equation. In this case, the wave equation can analytically be solved in the slow-roll approximation. The Bunch-Davies vacuum condition in the asymptotic past is imposed because the modes lie well inside the Hubble radius. Of course, if a general vacuum condition is imposed, there are additional features in the power spectrum. During inflation, quantum fluctuations are amplified and stretched, and then nearly frozen on super-Hubble scales. The single-field slow-roll inflation predicts a slightly red-tilted spectrum: and a consistency relation n T = −r/8 between the tensor spectral index n T and the tensorto-scalar ratio r. Here, n T and r are evaluated at the epoch when a given perturbation mode leaves the Hubble horizon, i.e. k = aH. We usually introduce the tensor-to-scalar ratio r ≡ A T /A S at a pivot scale k * , i.e. the amplitude of tensor perturbations A T with respect to that of scalar perturbations A S . The red-tilted spectrum means that the amplitude of tensor perturbations becomes small on small scales due to the fact that large-scale modes are earlier stretched across the Hubble horizon than small-scale modes as the energy density slowly decreases during inflation. From the wave equation we see that the evolution of h ij depends explicitly on a(τ ) and implicitly on inflation potential through FRW equations. Hence the power spectrum of GWs encodes useful information on the evolution of the scale factor. Moreover, the tensor-to-scalar ratio is related to the energy scale of inflation by V = 3π 2 M 4 pl A s r/2 = (1.88 × 10 16 GeV) 4 r/0.10. Here we have adopted the estimated value of the amplitude of scalar perturbations from the Planck 2015 data [11]. Different models predict different values of the tensor-to-scalar ratio. For example, the simplest chaotic inflation with a quartic potential [10] predicts a large value of r ≈ 0.26 while the R 2 inflation [12] predicts a small value of r ≈ 0.0033 to lowest order in slow-roll parameters if the number of e-folds N * = 60 is assumed. With the help of the scalar spectral index, the estimated value of r is robust to discriminate slow-roll inflationary models. In summary, the measurement of the power spectrum of GWs helps us to • test the vacuum initial condition, • detect the evolution of the scale factor, • determine the energy scale of inflation, • discriminate inflationary models.
If the primordial GWs are detected, the next important question to answer is what is the shape of the power spectrum of GWs and whether there are additional features in the power spectrum. In the slow-roll inflationary scenario, the shape of the tensor power spectrum is characterized by the tensor spectral index n T since the running of the spectral index is negligible to the lowest order in slow-roll parameters. More general shapes beyond slow-roll may be reconstructed by using a binning method of a cubic spline interpolation in a logarithmic wavenumber space [13][14][15]. Checking the consistency relation between n T and r provides a powerful test of the single-field slow-roll inflationary scenario. The violation of the consistency relation could in principle come from the following two aspects. The first is that the second and third terms on the left hand side of Eq. (2.2) are modified in general inflationary models, such as the k-inflation [16], Gauss-Bonnet inflation [17][18][19] and generalized G-inflation [20]. In the model of k-inflation, the consistency relation becomes n T = −r/8c S , where c S is the sound speed of scalar perturbations. In the Gauss-Bonnet inflationary model, the consistency relation is broken due to n T = −r/8 − δ 1 , where δ 1 is determined by the Gauss-Bonnet coupling term. In the model of generalized G-inflation, the tensor spectral index not only depends on the evolution of the scale factor but also the higher-order derivative terms, which admits a blue-tilted power spectrum of -5 -GWs. The second is that there exists a non-negligible source term on the right hand side of Eq. (2.2) during inflation. In this case, the wave equation is solved by the Green's function method. Possible sources of generating GWs include first-order scalar perturbations [21], perturbations of the extra field such as the curvaton [22] and spectator field [23,24], and particle production during inflation [25].

During Preheating
In the inflationary scenario, at the end of inflation, the inflaton field begins to oscillate around the minimum of its potential. Such coherent oscillations produce elementary particles and eventually reheat the Universe. This process is called reheating [26]. The coupling between the inflaton field and other fields is necessarily tiny ensuring that reheating proceeds slowly. Preheating provides a more rapidly efficient mechanism for extracting energy from the inflation field by parametric resonance [27]. Such a process is so rapid that the produced particles are not in thermal equilibrium. The preheating leads to large and time-dependent inhomogeneities of the stress tensor that source a stochastic background of GWs [28]. Unlike GWs produced during inflation, they are generated and remain in the Hubble horizon until now. Clearly, their wavelengths are smaller than the Hubble radius at the time of GW production. Therefore, the peak frequency of this type of stochastic GWs is typically of order more than 10 3 Hz. Detecting such high-frequency GWs is particularly challenging. For example, for the φ 4 and φ 2 chaotic inflationary models, lattice simulations [29] show that preheating can lead to GWs with frequencies of around 10 6 ∼ 10 8 Hz and peak power of Ω gw h 2 ≈ 10 −9 ∼ 10 −11 at present [30]. For hybrid inflation, GWs cover a larger range of frequencies. The peak wavelength depends essentially on the coupling constant [31]. See [32] for the most recent review on the GWs from the inflationary era and preheating era.

Observational Implications
Observational constraints on GWs produced during inflation mainly come from the B-mode polarization of the CMB anisotropies. Such GWs generate a quadrupolar anisotropy of momentum m = 2 in the intensity field of photons at the epoch of recombination while scalar perturbations generate only a quadrupolar anisotropy of momentum m = 0. Importantly, only the quadrupolar anisotropy of momentum m = 2 causes the B-mode polarization. Therefore, the measurement of the B-mode polarization allows us to probe GWs produced during inflation. The Planck 2015 data give the 95% confidence level upper bound for the tensor-to-scalar ratio r 0.002 < 0.10 at the pivot scale k * = 0.002 Mpc −1 [11], which is improved to r 0.002 < 0.08 by adding the BKP cross-correlation likelihood [33]. With the help of the scalar spectral index n S , slow-roll inflationary models are discriminated in the r − n S plane, as shown in Fig. 2. For example, the inflationary model with a quartic potential [10] is strongly disfavored by recent CMB data. Future ground-based and space-based CMB experiments will provide high precision measurements of the B-mode polarization. Actually, constraints on GWs produced during inflation from measurements of the CMB anisotropies are limited to a narrow scale range of 10 −4 ∼ 10 −1 Mpc −1 , which corresponds to the change of e-folds number ∆N ≈ 8. This implies that it is impossible to detect the global shape of the tensor spectrum by using only CMB data. Laser interferometer gravitational-wave experiments provide the possibility of measuring the stochastic background at small scales. During inflation, quantum fluctuations of h ij are amplified and stretched across the Hubble horizon, and then nearly frozen on super-Hubble scales. After reentering the Hubble horizon at the epoch of radiation or matter dominance, tensor perturbations begin to oscillate with the amplitude damped by a factor a −1 . Unfortunately, for the nearly scale-invariant tensor spectrum predicted by the single-field slow-roll inflationary models, the energy spectrum in the frequency range of 10 −10 ∼ 10 3 Hz at present becomes too weak to be directly detected by pulsar timing array (PTA) experiments and laser interferometer experiments [34]. If it is a blue-tilted tensor spectrum predicted by inflationary models beyond slow-roll or caused by the source term, the direct detection might be possible in the future [35].
For GWs produced during preheating, since their wavelengths are smaller than the Hubble radius at the time of GWs production, the corresponding peak frequencies are typically of order more than 10 3 Hz. It is particularly challenging to detect such highfrequency GWs by ground-based laser interferometer experiments.

The Gravitational Waves From Phase Transitions
It is believed that our observable Universe has experienced several PTs (PTs), during which an unknown high-degree symmetry is subsequently relaxed down to the broken EW symmetry that is well described by the SM of particle physics. If the PT from a high-temperature symmetric phase to a low-temperature symmetry-broken phase is of first order, the true vacuum bubbles would be nucleated within the false vacuum, leading to the expanding, colliding and merging bubbles that generate a stochastic background of GWs (see [36] for -7 -recent review and [37] for the discussions of GWs from cosmic strings and domain walls, which will not be discussed in this review). The primary motivations to study the GWs from PTs are two-folds: First, the EWPT of SM is cross-over according to the current measurements of Higgs mass. Therefore, any detections of GWs from PTs would necessarily provide us a unique probe beyond the SM, which cannot be directly probed by the particle colliders in a foreseeable future; second, the stochastic backgrounds of GWs also consist of the GWs from the inflation era and the reheating era and the other cosmological defeats from PTs like cosmic strings and domain walls. Therefore, it would help us to extract the GWs of astrophysical sources from the stochastic backgrounds of GWs. However, the topics of detections are not discussed in this review, which should merit another paper to discuss how to detect the stochastic background of GWs and how to distinguish the signals of PTs from those signals of reheating or other cosmological defeats.

Bubble Nucleation
In the seminal papers [38,39], the first semiclassical description of vacuum decay in flat spacetime was developed for a single scalar field without derivative interactions. The vacuum decay is implemented through barrier penetration from the unstable false vacuum to the stable true vacuum, of which the field configuration is captured in the so-called bounce equation. The bounce equation describes the true vacuum bubbles nucleated during barrier penetration in the surrounding false vacuum with the probability per unit time per unit volume being of form The coefficient B was worked out in [38] as the on-shell Euclidean action of the bounce solution and the coefficient A was worked out in [39] to properly account for the quantum corrections 1 . In the very special case where the potential barrier is larger than the difference of energy density between false and true vacuums, a thin-wall approximation was proposed in [38] to evaluate the bounce action in a closed form consisting both contributions from the bubble interior and the bubble wall. The insight from thin-wall approximation provides us a physical picture of bubble nucleation that the released energy from spreading true vacuum into false vacuum goes to the acceleration of the bubble wall until a critical bubble with zero energy is formed. The size of such critical bubble was determined by the stationary point of the Euclidean action, which gives a rough cancelation between the energy from bubble interior and the energy from bubble wall.
The seminal works [38,39] were later extended to the case of vacuum decay in curved spacetime in [40] and thermal decay in flat spacetime in [41]. Although the presence of gravity does not change the general picture, it does lower the probability of vacuum decay by nucleating larger bubbles. In the extreme case of almost degenerated vacuums, gravity can even stabilize the false vacuum by preventing us from making a zero-energy bubble in the false vacuum, which in the absence of gravity can be simply done by expanding the bubble interior to balance the surface tension of bubble wall. When the temperature turns on, instead of an O(4)-symmetric bubble in Euclidean spacetime, an O(3)-symmetric bubble solution that is periodic in time direction with period of inverse temperature T −1 is expected, and the size of such critical bubble is determined by the maximum point of total energy. The tunneling rate was estimated in [41,42] as is estimated at the bounce profile φ B (r) of equation-of-motion, The effective potential usually consists of tree-level potential, zero-temperature corrections and finite-temperature corrections including the daisy resummation. Recently, a new thermal resummation procedure was proposed in [43], which makes it possible to match theories to an EFT at finite temperature. Future works can be carried out along this direction. The naive strategy of solving above bounce equation is the shooting algorithm, which is illustrated in Fig. 3. When multiple fields are involved at play, four approaches were introduced in the literatures, the early trial [44], the undamping/damping algorithm [45], the wildly adopted path deformation algorithm [46], and the recently proposed semi-analytic perturbative approach [47]. The generalization of [38] in the case of non-zero cosmological constant allowed for both false and true vacuums was carried out in [48]. The vacuum decay is enhanced in the presence of gravity when the average of two vacuums is positive and suppressed when both vacuums are negative. The vacuum decay for a positive false vacuum but negative average of two vacuums is enhanced or suppressed when the gravitational effect becomes important or not. There is currently no satisfactory description of thermal decay in curved spacetime; however, we will clarify in the future work that the gravity correction only scales as κ/(T R 3 ), and it is not important because the characteristic size of bubble R is much larger than the Planck length √ κ ∼ 1/ √ G for the characteristic scale of PT below the Planck scale.

Bubble Expansion
After the bubble is nucleated within false vacuum, the bubble wall will rapidly expand until approaching the speed of light [38] or colliding with other bubble walls. However, the realistic background is actually a thermal plasma full of relativistic particles 2 , which  Figure 3. The schematic illustration of shooting algorithm for the bounce equation, which is equivalent to a particle moving in the inverse potential with Hubble friction. The exit point of field profile φ exit presents the field value at the center of nucleated bubble, which can be found by the shooting algorithm with errors and trials of under/over shooting until an exact shooting is achieved.
Here, the bubble profile is rescaled by the temperature M φ with respect to the bubble size rescaled by the mass scale M R of effective potential.
will interact with bubble walls by some form of frictions. The central problem of bubble expansion in plasma is to figure out the interconnections among following quantities: the -10 -bubble wall velocity v w , the friction η on the wall from plasma, the strength α of PT that measure the released vacuum energy density with respect to the radiation energy density, and the efficiency factors κ φ and κ v that measure the capability of transferring the liberated vacuum energy into the bubble wall expansion and the bulk fluid motion, respectively. These quantities serve as the interface between the theoretical models of particle physics and the simulations of bubble collisions. After the bubble is nucleated within thermal plasma, it starts stopping growing when the friction balances the net pressure from inside and outside of the bubble wall. If the pressure from outside of bubble wall is larger than the inside pressure while fluid velocity from outside of bubble wall is smaller than the inside velocity, the bubble wall behaves as deflagration, and the opposite definitions for detonation. Both deflagration and detonation are further classified as weak, Jouguet and strong types according to the fluid velocity from inside of bubble wall being smaller, equal and larger than the speed of sound. In the paper [54], it was proposed that with Jouguet condition, the velocity of bubble wall v w would be given by a simple formula [54] expressed in terms of the strength of PT α alone, so does the efficiency factor κ v . The formula of Jouguet detonation was extensively used in the literatures at the early stage of studies of GWs from PTs despite the fact [55] that the Jouguet detonation can be unrealistic in the cosmological setup of PTs. Since the work [55], several parallel explorations have been considered as follows.
• Beyond Chapman-Jouguet condition. Generalizing model-independent parametrization equation [55] of friction, [56] was able to explore the full range of bubble wall velocity for both deflagrations and detonations, where the analytic approximations were found for both non-relativistic and ultra-relativistic wall velocities. In the stateof-art work of [57], it gives a unified picture and user-friendly fitting formulas for the dynamical regimes of bubble expansion. A different but more accurate approach was adopted in [58] right after [56] with microscopic considerations on the particle contents of specific models, which was revisited and improved in the subsequent work [59]. Future works can be carried out along this direction.
• Criteria for runaway bubble walls. Apart from those stationary solutions of bubble wall with terminal velocity, there exists a runaway solution when the friction is too small to prevent the wall from approaching the speed of light. A simple criterion was found in [60] that the bubble walls will runaway if the effective potential in the true vacuum remains deeper than the false vacuum even after replacing the thermal potential by its second-order Taylor expansion term in the false vacuum, but not vice versa. This was reformulated in [57] simply by comparing the strength α of PTs with respect to a critical value α ∞ . Hence, combined considerations of hydrodynamics and microphysics of specific models were later extended in [61] and [62] to the runaway regime. Future works can be carried out along this direction.
• Reconciliations with baryogenesis. Large GWs signals from PTs require fast-moving bubble walls to collide with each other, while baryogenesis scenarios need slow-moving bubble wall to have an effective diffusion process. In [58], it also opens the possibility -11 -that can reconcile the baryogenesis with GWs from PTs in the presence of fermions with large Yukawa coupling and heavier stabilizing bosons. Later on, [63] pointed out that the relevant velocity responsible for the baryogenesis is the relative velocity between the wall and the plasma, which can be much smaller than the wall velocity when the strength of PTs becomes stronger. Therefore, it is possible to simultaneously generate large GWs from PTs without jeopardizing the effectiveness of baryogenesis.

Bubble Percolation
After bubble nucleation and bubble expansion, bubble percolation starts with colliding bubbles until PT is completed. When the initial size of nucleated bubble can be neglected, the duration of PT can be roughly estimated by the mean radius of bubbles at collisions, which is characterized by a single parameter β. Both β and the previous mentioned α are evaluated at nucleation temperature T * , which is defined as the temperature at which the number of generated bubble per unit time per Hubble volume is of order one. The old picture of bubble percolation consists of two sources, namely, the colliding bubble walls and the turbulent motions of bulk fluids along with their associated magnetic field. However, the new picture of bubble percolation is added with an extra source from the sound waves of bulk motion as the result of bubble collision, which exist long after the bubble percolation.
• Bubble collisions. It was first realized by Witten in [64] that the QCD PTs might leave us with detectable GWs from violent bubble collisions with their peak frequency characterized by the size of bubbles when they collided and with their peak amplitude estimated by the relative size of bubbles with respect to Hubble horizon at collision. Later on, this insight was generalized by Hogan in [65] to the case of EW PTs. In a series of papers [66][67][68][69], the preliminary simulations were first implemented for bubble collisions to capture the general features of GWs from PTs. It is found in [66,67] that remarkably the spectrum of GWs from simulating the collision of two vacuum bubbles in Minkowski space only depends upon the grossest features of the bubble collisions, namely, the strength α and duration β of PTs, which is similar to the result of simulation for hundreds of vacuum bubbles [68]. As we have shown in Fig. 4, the envelope approximation was proposed in [66][67][68] that the GWs are mainly generated from the uncollided envelopes of colliding bubble walls and any GWs from the overlap region can be neglected. The extension to the thermal bubble collisions was carried out later in the simulation [69], where the Jouguet detonation mode was explicitly used and the turbulent motion in fluid stirred by bubble collisions was appreciated in this case. The state-of-the-art results of the GWs spectrum from PTs due to colliding bubble walls were settled down in [70], which although disagreed with earlier study [71], they finally achieved consensus in [72]. Recently, under the thin-wall and envelope approximations in flat background, it was claimed in [73] that the GWs spectrum can be estimated analytically without need of simulations. Future works can be carried out along this direction from both theoretical and numerical points of view.
-12 -R true vacuum false vacuum thin shell envelope Figure 4. The schematic illustration of envelope approximation that the GWs are mainly generated from the uncollided envelopes of colliding bubble walls and any GWs from the overlap region can be neglected.
• Turbulent MHD. The possibility of generating GWs from turbulent motion of bulk fluid was mentioned long time ago in [64] as remnant of bubble collisions, which was first estimated in [69] as Kolmogorov spectrum under quadrupole approximation. Apart from above turbulent velocity field, the fully ionized plasma could also generate the turbulent magnetic field under turbulent motion, which itself is also a source of GWs. The early analysis [74,75] involving MHD exhibits three problems [37]: large-scale problem (addressed in [76,77]), time-evolution problem (persisted in [77] and addressed in [76]) and dispersion-relation problem (corrected in [76,77]). The state-of-art result of the GWs spectrum from PTs due to non-helical MHD turbulence was analytically established in [78], where it ignored the circularly polarized GWs [79][80][81][82][83][84][85] from helical turbulence due to the macroscopic parity violation. The numerical simulations of relativistic MHD turbulence are certainly needed to make further confirmations for the above analytic results in the future.
• Sound waves. The possibility of generating GWs from sound waves was originally pointed out in [65], although, it left forgotten for a long time until recent revelation in [86], where envelope approximation for modeling the colliding bubble walls is abandoned and the GWs should be dominated by the overlapping sound waves in the bulk fluid. The breakthrough findings of [86] were further quantitatively understood in the updated simulations [87,88] and theoretically modeled in [89], before which the possibility of having an inverse acoustic cascade was investigated in [90], suggesting the potentially very strong enhancement of the sound wave density at small wave number. Future works can be carried out along this direction for larger simulations with more bubbles.

Gravitational Waves
The GWs from the strong first-order PTs are guaranteed through above three processes: bubble nucleation, bubble expansion and bubble percolation. The bubble nucleation requires a potential barrier in order to tunnel from the false vacuum to the true vacuum. The bubble expansion requires fast enough moving bubble walls in order to have strong signals. The bubble percolation requires efficient collisions in order to dissipate released vacuum energy into bulk fluid motions. The fitting formulas of GWs spectrums from numerical simulations are well summarized in [36] and can be straightforwardly applied to the particle physics models, where the strength α and duration β evaluated at nucleation temperature T * along with the bubble wall velocity v w are obtained from the microscopic particle physics models, whereas the efficiency factors κ φ and κ v are approximated from the fitting formulas in [57]. There leaves one last discussion of which particle physics models can exhibit first-order PT so that GWs can be generated. The SM with phenomenological Higgs mechanism crosses over from the high-temperature phase to low-temperature phase [91] if the Higgs boson mass is larger than the W boson mass. To have GWs from first-order PTs, one has to go beyond SM (BSM) 3 . Here, we give an incomplete list of these models with explicit discussions on GWs, which will be revisited in the near future if we want to construct the reliable templates in order to extract the GW signals from the stochastic background.
• Higher dimensional operators. The simplest example is to add a cubic term, which is usually expected in the supersymmetric extensions of SM in order to make a strong first-order PTs. By adopting the polynomial fitting formulae [93] for the bounce action from a general quartic potential with a cubic term, a semi-analytic calculation of the GW signal from the EWPT was carried out in [94]. The other important example is the dimension-six operator [95,96], of which the GW signals from PTs were first analyzed in [97] and revisited in [98] with full scope of one-loop effective potential at finite temperature. In both examples, the bubble walls with detonations or run-away behavior have been considered in [99]. We produce in Fig. 5 the spectrum of GWs from PTs that originated from the dimension-six operator, which will be discussed in details in a future work [100].
• Additional scalar sectors. The simplest and extensively studied example is the gauge singlet scalar extension of SM [59,[101][102][103][104][105][106][107], which can be naturally fitted into previously mentioned cubic term [94] and can be tested at future colliders by measuring the triple Higgs coupling precisely 4 . The other important example is the charged scalar under SM gauge group, of which the simplest realization is the two-Higgs-doublet model (2HDM). The GW from 2HDM was preliminarily analyzed in [109] and further studied in [36] for the case of CP-conservation and recently revisited in [110] for the case of CP-violation.
• Supersymmetric extensions. The capability of detecting GWs from PTs in Minimal Supersymmetric Standard Model (MSSM) and Next-to-Minimal Supersymmetric Standard Model (NMSSM) were first estimated in [111] and further explored in [112]. Both the PTs from MSSM [59] and NMSSM [97] were thought to be not strong enough to produce significant GW signals; however, the parameter space with strong first-order PTs has been identified for a modified version of MSSM [113] and a general version of NMSSM [62].
• Hidden dark sectors. The cosmological implications concerning with possible production of GWs from hidden dark sector were first explored in [114], and later discussed for light GeV scalar [115], vector thermal dark matter [116], UV-conformal dark sector [117], SU (N ) dark sectors with n f flavors [118], dark U (1) gauge complex scalar singlet [119], hidden sector with run-away bubble walls [120], and PTs involving successive hidden gauge symmetry breaking [121], dark matter asymmetry [122] and two-step transition [123]. The GWs from first-order PTs could be a unique probe for these dark hidden sectors.

The Gravitational Waves From Binary Systems
GW sources can be divided into two categories, one is deterministic and the other is stochastic. For example, the primordial GW produced by the quantum fluctuations and by PTs during the early Universe belongs to the stochastic one. Due to the intrinsic random character, we cannot predict the waveform produced by the stochastic sources. The deterministic sources include two types. One type is predictable, while the other one is not. For example, the supernova is believed producing GWs. But the dynamics of supernova is too complicated to predict the related GW form. Another example is the foreground noise of LISA produced by white dwarf-white dwarf binaries in our milky way galaxy [134][135][136]. Due to the large number of such binaries, the combination of these GW signals makes the waveform unpredictable. In this section, we focus on the predictable GW sources. For such sources, we can construct theoretical model for gravitational waveform before GW detection experiment [137]. Then when the experiment data are ready, we can use matched filtering data analysis techniques to improve the GW detection sensitivity. And more, the theoretical model can be used to realize the astronomy detection of the GW sources through matched filtering [138,139]. Binary systems are among the most important and major sources for GW detection projects including PTA, space-based laser interferometer detectors such as eLISA, Taiji and Tianqin projects, ground-based laser interferometer detectors such as Advanced LIGO, Advanced Virgo, Kagra and others. The GW frequency of binary system around merger state can be characterized by the ISCO (most Inner Stable Circular Orbit) frequency: where M is the total mass of the binary system, and we have used geometric unit with c = G = 1. When M = M , f isco ≈ 2000Hz. At the same time, binary systems belong to predictable GW sources [139]. A typical example of the gravitational waveform and the related frequency are shown in Fig. 6. Note the time scale in the plot is different for inspiral stage and the merger/ringdown stage. The component of the binary system can be white dwarf, neutron star, black hole and even quark star. If gravitational theory beyond general relativity is true, other objects including axion star [140], gravastar [141] and other exotics might also be the component of binary system. Based on current understanding of stellar evolution, the final fate of a star may be white dwarf, neutron star or black hole which is determined by the mass of the star. Binary systems can be formed through capture process or many body interaction. Due to the gravitational slingshot effect involved in three-body interaction, binary systems can be -16 - . Typical example of the gravitational waveform and the corresponding frequency. In this example, the binary system is a binary black hole with equal mass. And the two black holes are spinless. The whole evolution process of a binary black hole system includes inspiral, merger and ringdown. But the boundary to distinguish these three stages is fuzzy. Compared to the inspiral stage, the merger and ringdown stages are much shorter. In this plot, we use different time scale for inspiral and merger/ringdown stage in order to make the merger/ringdown stage clearer.
formed through many body interaction. Because a binary system is dynamically stable, binary systems are expected to be quite common in our Universe.
The direct detection of GWs opens the GW astronomy. Through this new window to our Universe, GWs will not only bring us many new observations on kinds of objects in the Universe, but also give us many insights on fundamental physical theory. For example, black hole is a mystery from quantum gravity theory view. Especially, black hole presents us an information loss puzzle. Is it possible the black hole horizon be replaced with other objects such as firewall or anything else? Or any experiments present evidence that black hole horizon really exists? Unfortunately, one cannot get observational proof of the existence of a black hole event horizon based on EM waves [142]. In contrast, it is possible to use GW observation to give a direct proof of the existence of a black hole event horizon [143]. Currently, these promising aims can be most possibly achieved by predictable GW sources, especially binary systems. In order to realize these aims, we have to construct accurate enough gravitational waveform model for binary systems. There are three kinds of methods available to treat binary systems. They are PN approximation, numerical relativity and black hole perturbation method which are widely used for early inspiral stage, plunge and merger stage, and post merger stage, respectively.

Post-Newtonian Approximation
Based on quadrupole formula [144], we can estimate the GWs related to a binary system: where (r, θ, φ) is the position of the observer with respect to the binary, M , R and Ω are the total mass, the separation and the mutual orbiting frequency of the two components of the binary. Here h + and h × correspond to two polarization modes of GWs. These two modes are related to the dynamical metric form as with (t, r, θ, φ) corresponding to the transverse-traceless gauge [145]. Here we need point out two points. First one is that the planer wave approximation has been taken in the above equation which is reasonable because the field point is very far away from the source. The second is that above estimation about h + and h × is only leading order of the PN approximation which is reasonable because the weak GW source moves slowly. For most realistic GW sources, these approximations will break down. Even for cases in which the approximation is reasonable, above approximation cannot satisfy the need of GW detection [146,147]. For example, when the two components of a binary system are some far away, they move slowly. So, this kind of situation corresponds to weak and moving slowly GW source. Although the above approximation is reasonable, the accuracy is not good enough. Then, more PN order corrections are needed to improve the accuracy [148]. The theoretical framework of the PN approximation for the processing of binary systems consists of two parts: the dynamics of binary components and the GWform. In order to construct the GW model, we need to solve the binary dynamics, and then put the solution into the waveform theory part to get the explicit GW model. In general, the dynamical equations of the PN approximation are a highly nonlinear system of ordinary differential equations. In that case, the numerical method has to be employed to get the solution. GW models in which the GWform is expressed in time domain include TaylorT2, TaylorT4 and others [149]. Through stationary phase approximation, the post-Newtonian equations can be transformed into a frequency domain problem, and fortunately an approximate analytical expression can be obtained. Such a model is typically represented by the TaylorF2 model. Corrected by the binary black hole spin dynamics, especially the effect of the orbital precession, the TaylorF2 model is modified and improved by single and double precession model [150,151]. The TaylorT2 model is replaced by the X model when the elliptic orbit of the binary system is considered [152]. The theoretical model of the TaylorF2 model is extended to the elliptic orbit binary system including the post-circular model [153] and the enhanced post-circular model [154,155].

Numerical Relativity
Applying numerical methods to solve the Einstein equation is the topic of numerical relativity. Currently, the numerical relativity can only deal with the problem of GW modeling -18 -for the plunge and merger stage of binary systems. Because numerical relativity solves the Einstein equation without any approximation up to numerical error, this feature makes numerical relativity a universal tool to address kinds of GW sources. In Fig. 7, we show an example of the evolution process of a binary black hole system simulated by numerical relativity. But the nature of diffeomorphism invariance of the general relativity brings numerical calculation special difficulties. The study of numerical relativity began in the 60s of the last century [156], but the numerical instability made the code collapse after a few calculation steps. In 1990s, the construction of LIGO hardware strongly demanded the GW source modeling. So tens of universities and research institutions in USA and Europe jointly launched a Binary black hole Grand Challenge Project. However, the stability problem of numerical relativity has not been solved by that project. Out of one's expectation, the stability problem of numerical relativity was solved for the first time by Pretorius in 2005 [157]. In the following 2006, Baker group and Campanelli group also independently solved the instability problem [158,159]. Till now, there are more than ten numerical relativity groups around the world that have solved the stability problem, including the Princeton University, California Institute of Technology, University of Jena, Max Planck institute, the academy of mathematics and systems science, CAS [160] and others.
It is worth pointing out that within kinds of stable numerical relativity methods, which tips and treatments are necessary and/or sufficient for stable computation are still an open question. From the viewpoint of the theory of partial differential equations, one can only apply the hyperbolicity analysis to linearized Einstein equation [161]. In practice, however, the success of Pretorius in 2005 depends much on the formalism of the Einstein equations and the adaptive mesh refinement technique. Regarding the formalism, the so-called BSSN equation and the generalized harmonic coordinate equation have been widely used in numerical relativity. Adaptive mesh refinement is a very effective method to deal with multiscale problems. At present, the parallel adaptive mesh refinement codes developed especially for Einstein equations include the BAM code (developed by Bruegmann), AMSS-NCKU code (developed by Zhoujian Cao) [162], PAMR code (developed by Pretorius) and the Carpet code (developed by Schnetter).
Besides the stability problem, the major issues involved in numerical relativity are accuracy and efficiency. As for the formalism factor, more and more investigations show that Z4c and CCZ4 formalism is better than the BSSN formalism [163][164][165]. As an example of partial differential equation formalism of Einstein equations, we show Z4c formalism as following, which is proposed first time in 3D form in [164]: The unknown functions which need to be solved numerically are listed in the left-hand side of the above equations. For the meaning of the notations used in the above equations, we refer our readers to [164]. Spectral method, finite difference method and finite element method are three categories of numerical methods for partial differential equations. The vast majority of numerical relativity groups take the finite differential method. The SpEC code of California Institute of Technology uses spectral method. On the other hand, the application of finite element method into numerical relativity is still very few [166]. The spectral method has good computational efficiency due to its exponential convergence. However, the characteristics of its global data exchange limits its ability of parallel computing. The finite difference method working with domain decomposition algorithm can achieve good parallel computing efficiency. In order to deal with the multiscale characteristics of astrophysics, the refinement of the multilayer data structure is essential in the finite difference method. However, this method limits the number of cores for parallel computation by a single data layer, thus limits the strong parallel scalability. In contrast, the finite element method can combine the exponential convergence of the spectral method and the high parallel scalability of the difference method. Therefore, in principle, it is expected that the finite element calculation of the Einstein equations can achieve good strong parallel scalability. However, it is still an issue to construct the weak form of the Einstein equations, and to realize large-scale scientific computation [137,139]. Now the major challenges to numerical relativity placed by GW source modeling are the calculation of binary systems with mass ratio more than 100 to 1 and the well-converged simulation of binary systems including neutron stars [138].

Black Hole Perturbation
After the merger of the binary components, the system will ringdown and finally settle down to a Kerr black hole. If the two components are white dwarf and/or neutron star, the final remnant could be a stable neutron star. Here, we only concern with the case with Kerr black hole as the final product. Then the ringdown stage can be described by perturbation of a Kerr black hole. The black hole perturbation theory was pioneered by Regge, Wheeler and Zerilli [167,168] in which perturbation around a Schwarzschild black hole was considered. Teukolsky investigated the perturbation of a Kerr black hole on spacetime curvature level [169,170]. The former is called the Regge-Wheeler-Zerilli equation, and the latter is called the Teukolsky equation. The perturbation method used by Teukolsky is also valid to Schwarzschild black hole. For Schwarzschild black hole, Regge-Wheeler-Zerilli method and Teukolsky method is equivalent [171]. In the early 70s of the last century, some scholars have begun to use the Regge-Wheeler-Zerilli equation to study the associated GW for a test particle falling into a black hole [172][173][174]. Because the Teukolsky equation admits spin of a black hole, it is valid more extensively than the Regge-Wheeler-Zerilli equation. In recent years, there have been many works using the Teukolsky equation to calculate the GWs.
The Teukolsky equation can be solved by the PN approximation method (note that the Teukolsky equation is the target in question here, instead of Einstein equation itself like what described in the above subsection) or through numerical method. The PN ap-proximation method was mainly developed by Mano, Suzuki, Takasugi and others. These authors developed some hypergeometric functions and Coulomb wave functions to approximate the homogeneous solution of the Teukolsky equation [175,176]. Later on, Fujita and his coworkers applied PN approximation method to these analytical solutions to calculate the GWform and the related energy flux for extreme mass ratio binary systems. Currently, 22PN accuracy for Schwarzschild black hole and 11PN accuracy for Kerr black hole have been achieved [177,178]. However, such accuracy still do not yet satisfy the requirement of GW detection experiment [179].
The numerical methods to solve the Teukolsky equation can be divided into two categories: frequency domain method and time domain method. The frequency domain method divides the original Teukolsky equation into two ordinary differential equations through variable separating method and the Fourier expansion [169]. Regarding the radial equation, we can transform it to Sasaki-Nakamura equation before solving it numerically [180]. The Teukolsky equation can only provide the energy flux and the GWs. In order to provide the orbit of the test particle, Hughes and his coworkers applied adiabatic approximation method to the geodesic equation to treat the inspiral [181]. Now this kind of approximation has become an important method to treat binary systems with extreme mass ratio. Different from the usual finite difference numerical method, Fujita and Tagoshi used the hypergeometric function and the Coulomb wave function to expand the original equations which gives a quicker and more accurate numerical method [182]. The numerical method in time domain needs to solve a set of 2+1 partial differential equations. At present, one mainly uses the finite difference method to solve the problem (but see [183] for the finite element method). Due to the extreme mass ratio of the involved binary system, the related GW signal may last several years. Although the computation required by Teukolsky equation is much simpler than numerical relativity, how to improve the computational efficiency of the Teukolsky equation is also a great challenge for GW source modeling.

Cosmological Probes
In 1986, Schutz showed that from the GWs one can infer the Hubble constant by using the fact that the GWs from the binary systems encode the absolute distance information [184]. The inspiraling and merging compact binaries consisting of neutron stars and black holes can be considered as standard candles, or "standard sirens". The name of "siren" is due to the fact that the GW detectors are omni-directional and detect coherently the phase of the wave, which makes them in many ways more like microphones for sound than like conventional telescopes. For an expanding Universe, the chirp waveform of the GW can be generalized to the cosmological case by multiplying all masses by the factor 1 + z and the physical distance D can be replaced by the luminosity distance d L [185,186]. In fact, the waveform h produced by compact binary inspirals during the inspiral phase is theoretically well described by the analytical solution: 11) which is valid at the lowest (Newtonian) order for the "cross" GW polarization (the same expression holds for the "plus" polarization with a difference dependence on the orientation of the binary's orbital plane ι). Here, M c (z) is the (redshifted) chirp mass, f the GW frequency at the observer and Φ(f ) its phase. For the standard sirens, the luminosity distance d L (z) is the most important parameter entering the waveform (4.11). Parameter estimation over the observed GW signal (4.11) can directly yield the value of the luminosity distance of the GW source, together with an uncertainty due to the detector noise. Once the signal is detected and characterized, the luminosity distance of the source can be extracted. This indicates that one can measure the luminosity distance directly without the need of the cosmic distance ladder: standard sirens are "self-calibrating". However, it is not possible to infer the redshift z from the GW signal of the binary of black holes because all of the observed parameters such as the masses and distance are redshifted by the same factor 1 + z. To use the distance measurement for cosmography, one has to obtain the redshifts of the GW sources by some independent methods.
Before considering how to obtain the redshift information, one should think about how accurately the distance can be measured. The performance of the standard siren is limited by several effects. Firstly, the intrinsic measurement uncertainty in the amplitude of the detector's response is simply the inverse of the signal-to-noise ratio (SNR) [187,188], which is related to the detector sensitivity. Second, the weak gravitational lensing will distort the measurement of the luminosity distance, producing a magnification or demagnification on the order of a few percent [189][190][191][192]. Furthermore, the largest contribution of the uncertainty comes from the limited direction and source orientation sensitivity, and there is a large correlation between the distance, the sky position and orientation of the source. Thus, a network of detectors is needed to measure the position and orientation of the binary and to break the degeneracy among these parameters. By observing a simultaneous detection of a beamed EM signal, one can determine the sky position accurately and also will help to improve the measurements of the distance. The EM signal, called the EM counterpart, is also an important issue to the GW on cosmography.
The most traditional way to obtain the redshift of a GW event is through an accompanying EM signal, the EM counterpart. The binary merger of an NS with either an NS (BNS) or BH (BHNS) is hypothesized to be the progenitor of a short and intense burst of γ-rays, a so-called SGRB [193]. An EM counterpart like SGRB can provide the redshift information if the host galaxy of the event can be pinpointed. Moreover, SGRBs are likely to be strongly beamed phenomena [194], which allow one to constrain the inclination of the compact binary system, breaking the distance-inclination degeneracy. The GWs with short γ-rays bursts or other EM counterparts as the standard sirens have been studied in various papers (see [187,188,190,191,[195][196][197][198][199][200][201][202][203][204][205][206][207][208][209][210], and references therein). For example, Nissanke et.al. [202] used MCMC method and found that a network of advanced LIGO detectors can constrain the Hubble constant to a 5% accuracy. In Ref. [203], the authors demonstrate that with 1000 GW events detected by ground-based Einstein Telescope 5 it is possible to constrain the Hubble constant and dark matter energy parameter up to ∆h 0 ∼ 5 × 10 −3 5 Website for ET: www.et-gw.eu -23 -and ∆Ω m ∼ 0.02 using Fisher matrix approach, which was also found by Cai et al. [188] using MCMC method. Furthermore, Cai et al. used the Gaussian Process method and found that the equation of state of the dark energy can be constrained to ∆w(z) ∼ 0.03 in the low redshift region, which gives a better constraint than Ref. [203]. For the space-based detector LISA 6 , the expansion of the Universe and interacting dark energy have also been studied by Refs. [208,210].
There also exist other methods to infer the redshift information of the GW events, such as galaxy catalog proposed by Schutz [184], neutron star mass distribution [195,211], and the tidal deformation of neutron stars [212,213]. Measuring the redshift associated with a GW event is one of the biggest challenges in the future, see [214]; however, for a new GW standard sirens probe without redshift information by utilizing those BH binaries to be a tracer of the large-scale structure [215]. In addition, the spin of BH can also help us estimate the GW's parameters. GWs as the standard sirens to probe the cosmological parameters can provide an independent and complementary alternative to current experiments. It is expected that combining GW data with other astronomical observations such as supernovae, and adopting a better data analysis approach, the cosmological parameters could be constrained more precisely than the current situation.

Conclusions
The direct detection of GWs by LIGO initiates a new era of GW astronomy and GW cosmology. The GW physics is not only related to gravitational physics, but also closely related to particle physics, cosmology and astrophysics. The GWs provide us a new powerful tool to reveal various secrets of the nature. Indeed, a lot of relevant papers have appeared since the announcement of the direct detection of GWs. In this paper, we have briefly introduced three kinds of GW sources and relevant physics. They are GWs produced during inflation and preheating in the early Universe, from cosmic PTs and dynamics of compact binary systems, respectively. We also have discussed in a simple way the GWs as standard siren in the evolution of the Universe. Due to the limitation of space, we are not able to discuss all aspects of GW physics, but only focus on some main issues. Of course, it is also impossible to list a complete list of the references, quite probably some important references are missed here, we should apologize for this.