Note on gauge invariance of second order cosmological perturbations

  • We study the gauge invariant cosmological perturbations up to the second order. We demonstrate that there are infinite families of gauge invariant variables at both the first and second orders. The conversion formulae among different families are verified to be described by a finite number of bases that are gauge invariant. For the second order cosmological perturbations induced by the first order scalar perturbations, we explicitly represent their equations of motion in terms of the gauge invariant Newtonian, synchronous and hybrid variables, respectively.
  • 加载中
  • [1] J. M. Bardeen, Phys. Rev. D 22, 1882 (1980) doi: 10.1103/PhysRevD.22.1882
    [2] G. F. R. Ellis and M. Bruni, Phys. Rev. D 40, 1804 (1989)
    [3] M. Bruni, P. K. S. Dunsby, and G. F. R. Ellis, Astrophys. J. 395, 34 (1992) doi: 10.1086/171629
    [4] C. F. Sopuerta, M. Bruni, and L. Gualtieri, Phys. Rev. D 70, 064002 (2004), arXiv:gr-qc/0306027 doi: 10.1103/PhysRevD.70.064002
    [5] M. Bruni, S. Matarrese, S. Mollerach et al., Class. Quant. Grav. 14, 2585 (1997), arXiv:gr-qc/9609040[gr-qc doi: 10.1088/0264-9381/14/9/014
    [6] S. Sonego and M. Bruni, Commun. Math. Phys. 193, 209 (1998), arXiv:gr-qc/9708068 doi: 10.1007/s002200050325
    [7] K. Nakamura, Phys. Rev. D 74, 101301 (2006), arXiv:grqc/0605107[gr-qc doi: 10.1103/PhysRevD.74.101301
    [8] V. F. Mukhanov, H. Feldman, and R. H. Brandenberger, Phys. Rept. 215, 203 (1992) doi: 10.1016/0370-1573(92)90044-Z
    [9] P. A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020)
    [10] S. Matarrese, S. Mollerach, and M. Bruni, Phys. Rev. D 58, 043504 (1998), arXiv:astro-ph/9707278 doi: 10.1103/PhysRevD.58.043504
    [11] K. A. Malik and D. Wands, Phys. Rept. 475, 1 (2009), arXiv:0809.4944[astro-ph doi: 10.1016/j.physrep.2009.03.001
    [12] G. Domènech and M. Sasaki, Phys. Rev. D 97, 023521 (2018), arXiv:1709.09804[gr-qc doi: 10.1103/PhysRevD.97.023521
    [13] J.-O. Gong, (2019), arXiv: 1909.12708[gr-qc]
    [14] D. H. Lyth and Y. Rodriguez, Phys. Rev. D 71, 123508 (2005), arXiv:astro-ph/0502578 doi: 10.1103/PhysRevD.71.123508
    [15] M. Bruni, R. Crittenden, K. Koyama et al., Phys. Rev. D 85, 041301 (2012), arXiv:1106.3999[astro-ph.CO
    [16] M. Bruni and S. Sonego, Class. Quant. Grav. 16, L29 (1999), arXiv:gr-qc/9906017 doi: 10.1088/0264-9381/16/7/101
    [17] J. Yoo and R. Durrer, JCAP 09, 016 (2017), arXiv:1705.05839 [astro-ph.CO
    [18] K. Nakamura, Class. Quant. Grav. 31, 135013 (2014), arXiv:1403.1004[gr-qc doi: 10.1088/0264-9381/31/13/135013
    [19] S. Mollerach, D. Harari, and S. Matarrese, Phys. Rev. D 69, 063002 (2004), arXiv:astro-ph/0310711 doi: 10.1103/PhysRevD.69.063002
    [20] K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D 75, 123518 (2007), arXiv:gr-qc/0612013 doi: 10.1103/PhysRevD.75.123518
    [21] D. Baumann, P. J. Steinhardt, K. Takahashi et al., Phys. Rev. D 76, 084019 (2007), arXiv:hep-th/0703290 doi: 10.1103/PhysRevD.76.084019
    [22] H. Assadullahi and D. Wands, Phys. Rev. D 81, 023527 (2010), arXiv:0907.4073[astro-ph.CO doi: 10.1103/PhysRevD.81.023527
    [23] S. Saga, K. Ichiki, and N. Sugiyama, Phys. Rev. D 91, 024030 (2015), arXiv:1412.1081[astro-ph.CO doi: 10.1103/PhysRevD.91.024030
    [24] L. Alabidi, K. Kohri, M. Sasaki et al., JCAP 09, 017 (2012), arXiv:1203.4663[astro-ph.CO
    [25] L. Alabidi, K. Kohri, M. Sasaki et al., JCAP 05, 033 (2013), arXiv:1303.4519[astro-ph.CO
    [26] K. Inomata and T. Nakama, Phys. Rev. D 99, 043511 (2019), arXiv:1812.00674[astro-ph.CO doi: 10.1103/PhysRevD.99.043511
    [27] N. Orlofsky, A. Pierce, and J. D. Wells, Phys. Rev. D 95, 063518 (2017), arXiv:1612.05279[astro-ph.CO doi: 10.1103/PhysRevD.95.063518
    [28] I. Ben-Dayan, B. Keating, D. Leon et al., JCAP 06, 007 (2019), arXiv:1903.11843[astro-ph.CO
    [29] K. Kohri and T. Terada, Phys. Rev. D 97, 123532 (2018), arXiv:1804.08577[gr-qc doi: 10.1103/PhysRevD.97.123532
    [30] J. R. Espinosa, D. Racco, and A. Riotto, JCAP 09, 012 (2018), arXiv:1804.07732[hep-ph
    [31] S. Wang, T. Terada, and K. Kohri, Phys. Rev. D 99, 103531 (2019), arXiv:1903.05924[astro-ph.CO doi: 10.1103/PhysRevD.99.103531
    [32] S. J. Kapadia, K. L. Pandey, T. Suyama et al., Phys. Rev. D 101, 123535 (2020), arXiv:2005.05693[astro-ph.CO doi: 10.1103/PhysRevD.101.123535
    [33] S. J. Kapadia, K. L. Pandey, T. Suyama et al., (2020), arXiv: 2009.05514[gr-qc]
    [34] R.-g. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122, 201101 (2019), arXiv:1810.11000[astro-ph.CO doi: 10.1103/PhysRevLett.122.201101
    [35] C. Yuan, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D 101, 043019 (2020), arXiv:1910.09099[astro-ph.CO doi: 10.1103/PhysRevD.101.043019
    [36] C. Yuan and Q.-G. Huang, (2020), arXiv: 2007.10686[astroph.CO]
    [37] E. Bugaev and P. Klimai, Phys. Rev. D 81, 023517 (2010) doi: 10.1103/PhysRevD.81.023517
    [38] E. Bugaev and P. Klimai, Phys. Rev. D 83, 083521 (2011), arXiv:1012.4697[astro-ph.CO doi: 10.1103/PhysRevD.83.083521
    [39] R. Saito and J. Yokoyama, Phys. Rev. Lett. 102, 161101 (2009), [Erratum: Phys. Rev. Lett.107, 069901(2011)], arXiv: 0812.4339[astro-ph]
    [40] R. Saito and J. Yokoyama, Prog. Theor. Phys. 123, 867 (2010), [Erratum: Prog. Theor. Phys.126, 351(2011)], arXiv: 0912.5317[astro-ph.CO]
    [41] T. Nakama and T. Suyama, Phys. Rev. D 94, 043507 (2016), arXiv:1605.04482[gr-qc doi: 10.1103/PhysRevD.94.043507
    [42] Y.-F. Cai, C. Chen, X. Tong et al., Phys. Rev. D 100, 043518 (2019), arXiv:1902.08187[astro-ph.CO doi: 10.1103/PhysRevD.100.043518
    [43] K. Inomata and T. Terada, Phys. Rev. D 101, 023523 (2020), arXiv:1912.00785[gr-qc doi: 10.1103/PhysRevD.101.023523
    [44] C. Yuan, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D 101, 063018 (2020), arXiv:1912.00885[astro-ph.CO doi: 10.1103/PhysRevD.101.063018
    [45] M. Giovannini, (2020), arXiv: 2005.04962[hep-th]
    [46] J.-C. Hwang, D. Jeong, and H. Noh, Astrophys. J. 842, 46 (2017), arXiv:1704.03500[astro-ph.CO doi: 10.3847/1538-4357/aa74be
    [47] Y. Lu, A. Ali, Y. Gong et al., (2020), arXiv: 2006.03450[gr-qc]
    [48] V. De Luca, G. Franciolini, A. Kehagias et al., JCAP 2003, 014 (2020), arXiv:1911.09689[gr-qc
    [49] K. Tomikawa and T. Kobayashi, Phys. Rev. D 101, 083529 (2020), arXiv:1910.01880[gr-qc doi: 10.1103/PhysRevD.101.083529
    [50] E. Bertschinger, in Les Houches Summer School on Cosmology and Large Scale Structure (Session 60) (1993) pp. 273–348, arXiv: astro-ph/9503125
    [51] E. Lifshitz, J. Phys. (USSR) 10, 116 (1946)
    [52] K. Nakamura, PTEP 043, (2013), arXiv:1105.4007[gr-qc
    [53] K. Nakamura, Prog. Theor. Phys. 117, 17 (2007), arXiv:grqc/0605108 doi: 10.1143/PTP.117.17
    [54] K. Nakamura, (2020), 10.9734/bpi/taps/v3, arXiv: 1912.12805 [gr-qc]
    [55] M. Hohmann, Phys. Rev. D 101, 024061 (2020), arXiv:1910.09245[gr-qc doi: 10.1103/PhysRevD.101.024061
    [56] C. Uggla and J. Wainwright, Class. Quant. Grav. 29, 105002 (2012), arXiv:1112.0880[gr-qc doi: 10.1088/0264-9381/29/10/105002
    [57] K. A. Malik and D. Wands, (1998), arXiv: gr-qc/9804046
    [58] K. Inomata, (2020), arXiv: 2008.12300[gr-qc]
    [59] C. Pitrou, X. Roy, and O. Umeh, Class. Quant. Grav. 30, 165002 (2013), arXiv:1302.6174[astro-ph.CO doi: 10.1088/0264-9381/30/16/165002
    [60] S. Weinberg, Cosmology (2008)
    [61] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology (Oxford University Press, 2018)
    [62] J. W. York, Jr., Annales de L’Institut Henri Poincare Section (A) Physique Theorique 21, 319 (1974)
  • 加载中

Get Citation
Zhe Chang, Sai Wang and Qing-Hua Zhu. Note on gauge invariance of second order cosmological perturbations[J]. Chinese Physics C. doi: 10.1088/1674-1137/ac0c74
Zhe Chang, Sai Wang and Qing-Hua Zhu. Note on gauge invariance of second order cosmological perturbations[J]. Chinese Physics C.  doi: 10.1088/1674-1137/ac0c74 shu
Milestone
Received: 2020-12-14
Article Metric

Article Views(1574)
PDF Downloads(43)
Cited by(0)
Policy on re-use
To reuse of subscription content published by CPC, the users need to request permission from CPC, unless the content was published under an Open Access license which automatically permits that type of reuse.
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Email This Article

Title:
Email:

Note on gauge invariance of second order cosmological perturbations

Abstract: We study the gauge invariant cosmological perturbations up to the second order. We demonstrate that there are infinite families of gauge invariant variables at both the first and second orders. The conversion formulae among different families are verified to be described by a finite number of bases that are gauge invariant. For the second order cosmological perturbations induced by the first order scalar perturbations, we explicitly represent their equations of motion in terms of the gauge invariant Newtonian, synchronous and hybrid variables, respectively.

    HTML

    I.   INTRODUCTION
    • The gauge invariance is known as an important concept in the cosmological perturbation theory. In the first order, the gauge invariant variables were introduced for the first time by Bardeen in 1980s [1]. Since then, different gauge-invariant formalisms have been studied [2-7]. In addition, they have been widely utilized to study cosmological fluctuations [8], for example, the anisotropies in the cosmic microwave background and the large scale structures, which have been precisely measured in the past decades [9]. For the second order cosmological perturbations, gauge invariance has also been recently studied, and a quantity of gauge invariant variables have been constructed [5-7, 10-17]. In addition, it is believed that the gauge invariance is available to study higher order cosmological perturbations [18].

      However, recent changes were noticed when the second order cosmological perturbations were studied, including but not limited to the second order gravitational waves [19-23], which are induced by the first order scalar perturbations [24-42]. The energy density spectrum of the second order gravitational waves, as a physical observable, have recently been involved in the gauge issue [43-49]. A possible explanation is related to gauge fixing, which may trigger unknown fictitious perturbations [45]. Therefore, this problem can be addressed by constructing the gauge invariant variables [2-7, 10-13]. The other explanation is related to the definition of the physical observable, which has been suggested to be defined in the synchronous frame [48]. However, in such a framework, the concept of gauge invariance was suggested to be necessarily abandoned, because it is impossible to truly construct the gauge invariant second order synchronous variables [50].

      However, we would demonstrate that the gauge invariance can be preserved when second order cosmological perturbations are studied, in particular, the second order gravitational waves. On one hand, the power spectrum as the physical observable should be gauge invariant. Otherwise, it could take an arbitrary value if we choose an appropriate gauge fixing, e.g., the synchronous gauge [51], as reviewed in Ref. [50]. On the other hand, the gauge invariant synchronous variables can be logically defined for the induced gravitational waves at second order, though they are poorly defined for the second order scalar and vector perturbations. Accordingly, the energy density spectrum of second order gravitational waves could be well defined.

      In this study, we investigate the gauge invariance of second order cosmological perturbations by following the Lie derivative method [5-7, 10, 11]. The gauge invariant variables of any order can be systematically constructed in such a method. We will study the gauge invariant Newtonian (i.e., Bardeen's) and synchronous variables at first and second orders. In particular, for the first time, we will construct the so-called gauge invariant hybrid variables, i.e., Newtonian at first order while synchronous at second order, and vice versa. Further, we will determine the conversion formulae among different families of gauge invariant variables. Finally, we will derive the gauge invariant equations of motion for the second order cosmological perturbations by adopting the scalar-vector-tensor decomposition [52]. To accomplish this objective, we will decompose the gauge invariant perturbed Einstein field equations into the scalar, vector, and tensor components.

      The remainder of the paper is arranged as follows. In Section II, we introduce the gauge transformations and the gauge invariant variables, in accordance with the Lie derivative method. In Section III, we revisit the gauge invariant first order variables in such a method. In Section IV, we present explicit expressions of the gauge invariant second order variables, and study the conversion formulae among different families of gauge invariant variables. In Section V, the equations of motion for the cosmological perturbations are presented in terms of the gauge invariant Newtonian, synchronous and hybrid variables. Finally, the conclusions and discussions are summarized in Section VI.

    II.   GAUGE TRANSFORMATION AND GAUGE INVARIANT VARIABLES

      A.   Gauge transformation of an arbitrary tensor

    • A gauge transformation for a perturbed quantity in space-time starts from an infinitesimal transformation of coordinate $ x^{\mu} \rightarrow \tilde{x}^{\mu} $. The expansion of $ \tilde{x}^{\mu} $ up to second order is given by

      $ \tilde{x}^{\mu} = x^{\mu} + \xi^{(1) \mu} + \frac{1}{2} \xi^{(2) \mu} +{\cal{O}} \left(\xi^{(3)}\right) ,$

      (1)

      where $ \xi^{(1) \mu} $ and $ \xi^{(2) \mu} $ are the first and second order expansions of $ \tilde{x}^{\mu} $, respectively. For a generic tensor $ {\cal{Q}} $, to second order, the infinitesimal transformations are shown to be [5]

      $ \tilde{{\cal{Q}}}^{(1)} = {\cal{Q}}^{(1)} +{\cal{L}}_{\xi_1} {\cal{Q}}^{(0)}, \tag{2a}$

      $ \tilde{{\cal{Q}}}^{(2)} = {\cal{Q}}^{(2)} + 2{\cal{L}}_{\xi_1} {\cal{Q}}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) {\cal{Q}}^{(0)}, \tag{2b}$

      where we let $ \xi_1^{\mu} \equiv \xi^{(1) \mu} $ and $ \xi_2^{\mu} \equiv \xi^{(2) \mu} - \xi^{(1)\nu} \partial_{\nu} \xi^{(1) \mu} $, $ {\cal{Q}}^{(n)} $ denotes the n-th order perturbation of the tensor $ {\cal{Q}} $, and $ {\cal{L}}_{\xi} $ is Lie derivative along the infinitesimal vector $ \xi $. The degree of freedom of the transformation is the same for the first and second order perturbations. The $ \xi_1 $ determines $ \tilde{{\cal{Q}}}^{(1)} $. Once $ \xi_1 $ is fixed, $ \xi_2 $ determines the $ \tilde{{\cal{Q}}}^{(2)} $ . These formulae are based on the fact that the $ \xi_1 $ i is set to be in the same order of tensor perturbation $ {\cal{Q}}^{(1)} $. It could simplify the formalism of the perturbation theory. For a general infinitesimal transformation, there is no relevance between $ \xi_1 $ and $ {\cal{Q}}^{(1)} $ on the orders. In Appendix A, an introduction to the Lie derivative is shown. We briefly summarize the derivations of Eqs. (2) in Appendix B.

    • B.   Gauge invariant metric perturbations

    • Based on Eqs. (2a) and (2b), the infinitesimal transformations of the first and second order metric perturbations are presented as

      $ \tilde{g}_{\mu \nu}^{(1)} = g_{\mu \nu}^{(1)} +{\cal{L}}_{\xi_1} g_{\mu \nu}^{(0)},\tag{3a}$

      $ \tilde{g}^{(2)}_{\mu \nu} = g_{\mu \nu}^{(2)} + 2{\cal{L}}_{\xi_1} g_{\mu \nu}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) g_{\mu \nu}^{(0)} . \tag{3b}$

      If we introduce the gauge invariant metric perturbations via $ g_{\mu \nu}^{(\mathrm{GI}, i)} \equiv g_{\mu \nu}^{(i)} - \mathbb{C}_{\mu \nu}^{(i)} $ ($ i = 1, 2 $a for example), by adopting Eqs. (3a) and (3b), we can rewrite the counter terms $ \mathbb{C}_{\mu \nu}^{(i)} $ in terms of the infinitesimal vectors $ X^{\mu} $ and $ Y^{\mu} $, namely,

      $ g_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} = g_{\mu \nu}^{(1)} -{\cal{L}}_X g_{\mu \nu}^{(0)}, \tag{4a}$

      $ g^{(\mathrm{\rm{GI}, 2)}}_{\mu \nu} = g_{\mu \nu}^{(2)} - 2{\cal{L}}_X g_{\mu \nu}^{(1)} - ({\cal{L}}_Y -{\cal{L}}_X^2) g_{\mu \nu}^{(0)},\tag{4b} $

      where

      $ \tilde{X}^{\mu} = X^{\mu}+ \xi_1^{\mu}, \tag{5a}$

      $ \tilde{Y}^{\mu} = Y^{\mu}+ \xi_2^{\mu} + [\xi_1, X]^{\mu} . \tag{5b}$

      We present the derivations of Eqs. (4a)–(5b) in Appendix C. These formulae can also be found in Refs. [7, 53, 54]. Recently, they have been used in cosmology [7, 54] and in the post-Newtonian formalism [55].

      Based on Eqs. (5a) and (5b), the infinitesimal vectors $ X^{\mu} $ and $ Y^{\mu} $ could be independent of the metric perturbations in principle. In the first order, Bardeen constructed the gauge invariant variables in terms of the metric perturbations [1]. Therefore, we limit our investigations to the case in which both $ X^\mu $ and $ Y^\mu $ are expressed in terms of the metric perturbations in the following.

    • C.   Gauge invariant perturbations of the energy-momentum tensor

    • In the aspect of matter perturbations, we adopt the infinitesimal transformation in Eqs. (2a) and (2b) to obtain the perturbed energy-momentum tensors, i.e.,

      $\tilde{T}_{\mu \nu}^{(1)} = T_{\mu \nu}^{(1)} +{\cal{L}}_{\xi_1} T_{\mu \nu}^{(0)}, \tag{6a}$

      $ \tilde{T}^{(2)}_{\mu \nu} = T_{\mu \nu}^{(2)} + 2{\cal{L}}_{\xi_1} T_{\mu \nu}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) T_{\mu \nu}^{(0)} .\tag{6b} $

      The gauge invariant perturbed energy-momentum tensors can be given by

      $ T_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} = T_{\mu \nu}^{(1)} -{\cal{L}}_X T_{\mu \nu}^{(0)}, \tag{7a}$

      $ T^{(\mathrm{\rm{GI}, 2)}}_{\mu \nu} = T_{\mu \nu}^{(2)} - 2{\cal{L}}_X T_{\mu \nu}^{(1)} - ({\cal{L}}_Y -{\cal{L}}_X^2) T_{\mu \nu}^{(0)}, \tag{7b}$

      In the above formulae, it can be observed that the energy-momentum tensors are not in a special status. In general, any of the gauge invariant tensors could be formulated in the same form as Eqs. (7a) and (7b), i.e.,

      $ {\cal{Q}}^{(\mathrm{\rm{GI}, 1)}} = {\cal{Q}}^{(1)} -{\cal{L}}_X {\cal{Q}}^{(0)}, \tag{8a}$

      $ {\cal{Q}}^{(\mathrm{\rm{GI}, 2)}} = {\cal{Q}}^{(2)} - 2{\cal{L}}_X {\cal{Q}}^{(1)} - ({\cal{L}}_Y -{\cal{L}}_X^2) {\cal{Q}}^{(0)} . \tag{8b}$

    • D.   Gauge invariant Einstein field equations

    • It is already known that the gauge invariant perturbed Einstein field equations have the same form as the conventional ones [8, 54]. This can be explicitly expressed by expanding the Einstein field equations, $ G_{\mu \nu} = \kappa T_{\mu \nu} $, namely,

      $ \begin{aligned}[b] 0 =& G_{\mu \nu} - \kappa T_{\mu \nu} \\ \approx &G^{(0)}_{\mu \nu} - \kappa T_{\mu \nu}^{(0)} + G^{(1)}_{\mu \nu} - \kappa T_{\mu \nu}^{(1)} + \frac{1}{2} \left(G_{\mu \nu}^{(2)} - \kappa T_{\mu \nu}^{(2)}\right) \\ = &G_{\mu \nu}^{(0)} - \kappa T_{\mu \nu}^{(0)} + \left(G^{(\mathrm{\rm{GI}, 1)}} - \kappa T_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} +{\cal{L}}_X \left(G_{\mu \nu}^{(0)} - \kappa T_{\mu \nu}^{(0)}\right)\right) \\ &+ \frac{1}{2} \left(G_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} - \kappa T_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} + 2{\cal{L}}_X \left(G_{\mu \nu}^{(1)} - \kappa T_{\mu \nu}^{(1)}\right)\right. \\ &\left.+ \left({\cal{L}}_Y -{\cal{L}}_X^2\right) \left(G_{\mu \nu}^{(0)} - \kappa T^{(0)}_{\mu \nu}\right)\right), \end{aligned} $

      (9)

      where $ G_{\mu \nu} $ is the Einstein tensor, $ \kappa \equiv 8 \pi G $, and G is the gravitational constant. We can separate Eq. (9) in successive order to obtain

      $ G_{\mu \nu}^{(0)} = \kappa T_{\mu \nu}^{(0)}, \tag{10a}$

      $ G_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} = \kappa T_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}}, \tag{10b}$

      $ G_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} = \kappa T_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} . \tag{10c}$

      The gauge invariant perturbed Einstein tensor $ G_{\mu \nu}^{(\mathrm{\rm{GI}}, n)} $ takes the same form as the conventional one $ G_{\mu \nu}^{(n)} $. It can be obtained by substituting the metric perturbations $ g_{\mu \nu}^{(n)} $ with the gauge invariant ones $ g_{\mu \nu}^{(\mathrm{\rm{GI}}, n)} $. The same situation is also true for the perturbed energy-momentum tensor $ T_{\mu\nu}^{(\mathrm{\rm{GI}}, n)} $.

      In fact, the above conclusion is obvious, and can be understood as follows. As known, the perturbed Einstein tensors $ G_{\mu \nu}^{(n)} $ are defined with the metric perturbations $ g_{\mu\nu}^{(n)} $. Upon the gauge transformation, the transformed perturbed Einstein tensor $ \tilde{G}^{(n)}_{\mu \nu} $ should be written in terms of the transformed metric perturbations $ \tilde{g}_{\mu \nu}^{(n)} $. Based on this, because we notice that $ g_{\mu \nu}^{(\mathrm{\rm{GI}}, n)} $ in Eqs. (4a) and (4b) formally take the same form as the gauge transformations, i.e., $ \tilde{g}_{\mu \nu}^{(n)} $ in Eqs. (3a) and (3b), the gauge invariant perturbed Einstein tensor $ G_{\mu \nu}^{(\mathrm{\rm{GI}}, n)} $ can be defined with the gauge invariant metric perturbations $ g_{\mu\nu}^{(\mathrm{\rm{GI}}, n)} $.

    III.   GAUGE INVARIANT VARIABLES FOR THE FIRST ORDER COSMOLOGICAL METRIC PERTURBATIONS
    • The gauge invariant first order variables were first proposed by Bardeen [1]. As was known at the time, Bardeen's gauge invariant variables have the same form as the metric perturbations in the Newtonian gauge. In this section, we will reproduce Bardeen's formulae and further demonstrate that there are infinite families of gauge invariant variables that are allowed.

    • A.   Gauge transformations of the first order metric perturbations

    • In flat Friedmann-Lemaître-Robertson-Walker space-time, the metric takes the form of

      $ g_{\mu \nu}^{(0)} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} = a^2 (\eta) (- \mathrm{d} \eta^2 + \delta_{ij} \mathrm{d} x^i \mathrm{d} x^j),$

      (11)

      where $ \eta $ and $ a(\eta) $ are the conformal time and the scale factor of the Universe, respectively. The spatial curvature of the space-time is zero. The metric perturbations of the n-th order can take the form

      $ \begin{array}{l} g_{\mu \nu}^{(n)} = a^2 \left( \begin{array}{ccc} - 2 \phi^{(n)} & & \partial_i b^{(n)} + \nu_i^{(n)}\\ \partial_i b^{(n)} + \nu_i^{(n)} & & - 2 \psi^{(n)} \delta_{ij} + 2 \partial_i \partial_j {e}^{(n)} + \partial_i c_j^{(n)} + \partial_j c_i^{(n)} + h_{ij}^{(n)} \end{array} \right) , \end{array} $

      (12)

      where $ \phi^{(n)} $, $ \psi^{(n)} $, $ b^{(n)} $, $ {e}^{(n)} $ are scalar perturbations, $ \nu_i^{(n)} $ and $ c_j^{(n)} $ are vector perturbations, and $ h_{ij}^{(n)} $ are tensor perturbations. For tensor and vector perturbations, the transverse or traceless conditions should be satisfied as follows

      $ \partial_i \nu^{i (n)} = 0,\tag{13a}$

      $ \partial_i c^{i (n)} = 0, \tag{13b}$

      $ \delta^{jk} \partial_k h_{ij}^{(n)} = 0, \tag{13c}$

      $ \delta^{ij} h_{ij}^{(n)} = 0. \tag{13d}$

      These variables can be introduced via a scalar-vector-tensor decomposition, which is summarized in Appendix D.

      Based on Eq. (3a), the gauge transformations of the scalar, vector, and tensor perturbations are explicitly expressed as

      $ \begin{aligned}[b] \tilde{g}^{(1)}_{00} &= - 2 a^2 \tilde{\phi}^{(1)} \\ &= - 2 a^2 \phi^{(1)} - 2 a^2 \left( \partial_0 + \frac{\dot{a}}{a} \right) \xi_1^0,\end{aligned} \tag{14a}$

      $ \begin{aligned}[b] \tilde{g}^{(1)}_{0i} & = a^2 (\partial_i \tilde{b}^{(1)} + \tilde{\nu}_i^{(1)}) \\ & = a^2 (\partial_i b^{(1)} + \nu_i^{(1)}) + a^2 (\delta_{ij} \partial_0 \xi_1^j - \partial_i \xi^0_1), \end{aligned}\tag{14b} $

      $ \begin{aligned}[b] \tilde{g}^{(1)}_{ij} =& a^2 (- 2 \tilde{\psi}^{(1)} \delta_{ij} + 2 \partial_i \partial_j \tilde{e}^{(1)} + \partial_i \tilde{c}_j^{(1)} + \partial_j \tilde{c}_i^{(1)} + \tilde{h}_{ij}^{(1)}) \\ =& a^2 (- 2 \psi^{(1)} \delta_{ij} + 2 \partial_i \partial_j {e}^{(1)} + \partial_i c_j^{(1)} + \partial_j c_i^{(1)} + h_{ij}^{(1)}) \\ & + a^2 \left( (\delta_{ik} \partial_j + \delta_{jk} \partial_i) (\xi^k_{1, T} + \delta^{kl} \partial_l \xi_{1,S}) + \frac{2 \dot{a}}{a} \delta_{ij} \xi^0_1 \right) . \end{aligned}\tag{14c}$

      Here, the spatial part of $ \xi_1^{\mu} $ has been decomposed as $ \xi^i_1 = \xi^i_{1, T} + \delta^{ij} \partial_j \xi_{1,S} $, where we have the transverse $ \xi_{1, T}^i $ and longitudinal $ \xi_{1,S} $ parts. Using Eqs. (14a)-(14c) and the scalar-vector-tensor decomposition, we rewrite the gauge transformations of the metric perturbation variables as

      $ \tilde{\phi}^{(1)} = \phi^{(1)} + \left( \partial_0 + \frac{\dot{a}}{a} \right) \xi_1^0, \tag{15a}$

      $ \tilde{b}^{(1)} = b^{(1)} + \partial_0 \xi_{1,S} - \xi_1^0,\tag{15b} $

      $ \tilde{\nu}_i^{(1)} = \nu_i^{(1)} + \delta_{ij} \partial_0 \xi^j_{1, T}, \tag{15c}$

      $ \tilde{\psi}^{(1)} = \psi^{(1)} - \frac{\dot{a}}{a} \xi_1^0, \tag{15d}$

      $ \tilde{e}^{(1)} = {e}^{(1)} + \xi_{1,S}, \tag{15e}$

      $ \tilde{c}_i^{(1)} = c_i^{(1)} + \delta_{ik} \xi^k_{1, T}, \tag{15f}$

      $ \tilde{h}_{ij}^{(1)} = h_{ij}^{(1)} .\tag{15g} $

      In the following, we will introduce the gauge invariant variable of metric perturbations by adopting Eq. (15). Because the gauge invariant metric perturbations could take the same form as the Newtonian (or synchronous) gauge, we call them the gauge invariant Newtonian (or synchronous) metric perturbations for the purpose of presentations.

    • B.   Gauge invariant first order Newtonian variables

    • Based on Eq. (4a), we could obtain the gauge invariant Newtonian metric perturbations, i.e.,

      $ g_{00}^{(\mathrm{\rm{GI}, 1)}} = - 2 a^2 \Phi^{(1)}, \tag{16a}$

      $ g_{0 i}^{(\mathrm{\rm{GI}, 1)}} = a^2 V_i^{(1)}, \tag{16b}$

      $ g_{ij}^{(\mathrm{\rm{GI}, 1)}} = a^2 \left(- 2 \Psi^{(1)} \delta_{ij} + H_{ij}^{(1)}\right), \tag{16c}$

      where the gauge invariant variables are defined as

      $ \Phi^{(1)} = \phi^{(1)} - \left( \partial_0 + \frac{\dot{a}}{a} \right) X^0, \tag{17a}$

      $ \Psi^{(1)} = \psi^{(1)} + \frac{\dot{a}}{a} X^0, \tag{17b}$

      $ V^{(1)}_i = \nu_i^{(1)} - \delta_{ij} \partial_0 X^j_T, \tag{17c}$

      $ H_{ij}^{(1)} = h_{ij}^{(1)}, \tag{17d}$

      and we have decomposed $ X^i = : X_T^i + \delta^{ij} \partial_j X_S $, i.e., into the transverse and longitudinal parts. As expected, $ X^{\mu} $ is expressed in terms of the metric perturbations $ g^{(1)}_{\mu \nu} $, we can derive its formula from Eq. (16b),

      $ \begin{aligned}[b] a^2 V_i^{(1)} =& g_{0 i}^{(1)} -{\cal{L}}_X g_{0 i}^{(0)} \\ =& a^2 \left(\partial_i \left(b^{(1)} - \partial_0 X_S + X^0\right)\right. \\ &\left.+ \left(\nu_i^{(1)} - \delta_{ji} \partial_0 X^j_T\right)\right) . \end{aligned} $

      (18)

      Because the vector perturbation is transverse, it leads to

      $ \begin{array}{l} b^{(1)} - \partial_0 X_S + X^0 = 0. \end{array} $

      (19)

      By making use of Eq. (16c), namely,

      $ \begin{aligned}[b] a^2 (- 2 \Psi^{(1)} \delta_{ij} \!+\! H_{ij}^{(1)}) \!=& g_{ij}^{(1)} -{\cal{L}}_X g_{ij}^{(0)} \\ = &a^2 \Biggr(\! -\! 2 \delta_{ij} \left( \psi^{(1)} \!+\! \frac{\dot{a}}{a} X^0 \right) + 2 \partial_i \partial_j \left({e}^{(1)} \!-\! X_S\right) \\ &+ \left(\delta^k_j \partial_i + \delta^k_i \partial_j\right) \left(c_k^{(1)} - \delta_{lk} X^l_T\right) + h_{ij}^{(1)}\Biggr), \end{aligned} $

      (20)

      and the scalar-vector-tensor decomposition, we obtain

      $ \begin{array}{l} {e}^{(1)} - X_S = 0, \end{array} \tag{21a} $

      $ c_k^{(1)} - \delta_{lk} X_T^l = 0. \tag{21b}$

      Based on Eqs. (19), (21a) and (21b), $ X^{\mu} $ can be given by

      $ X^0 = \partial_0 {e}^{(1)} - b^{(1)}, \tag{22a} $

      $ X^i = \delta^{ik} \left(c_k^{(1)} + \partial_k {e}^{(1)}\right) . \tag{22b} $

      Therefore, we rewrite the gauge invariant variables as

      $ \Phi^{(1)} = \phi^{(1)} - \frac{1}{a} \partial_0 \left(a \left(\partial_0 {e}^{(1)} - b^{(1)}\right)\right), \tag{23a} $

      $ \Psi^{(1)} = \psi^{(1)} - \frac{\dot{a}}{a} \left(\partial_0 {e}^{(1)} - b^{(1)}\right), \tag{23b} $

      $ V_i^{(1)} = \nu_i^{(1)} - \partial_0 c_i^{(1)}, \tag{23c} $

      $ H_{ij}^{(1)} = h_{ij}^{(1)} . \tag{23d} $

      This implies that the gauge invariant variables proposed by Bardeen [1] can be reproduced by adopting the Lie derivative method [54].

    • C.   Gauge invariant first order synchronous variables

    • Based on Eq. (4a), the gauge invariant synchronous metric perturbations take the form of

      $ g_{00}^{(\mathrm{\rm{GI}, 1)}} = 0,\tag{24a} $

      $ g_{0 i}^{(\mathrm{\rm{GI}, 1)}} = 0, \tag{24b}$

      $ \begin{aligned}[b] g_{ij}^{(\mathrm{\rm{GI}, 1)}} =& a^2 \left(- 2 \Psi^{(1)} \delta_{ij} + 2 \partial_i \partial_j E^{(1)}\right.\\ &\left.+ \partial_i C_j^{(1)} + \partial_j C_i^{(1)} + H_{ij}^{(1)}\right), \end{aligned}\tag{24c} $

      where the gauge invariant variables are defined as

      $ \Psi^{(1)} = \psi^{(1)} + \frac{\dot{a}}{a} X^0, \tag{25a} $

      $ E^{(1)} = {e}^{(1)} - X_S, \tag{25b}$

      $ C_i^{(1)} = c_i^{(1)} - \delta_{ik} X_T^k, \tag{25c}$

      $ H_{ij}^{(1)} = h_{ij}^{(1)} . \tag{25d} $

      In this case, $ X^{\mu} $ can be determined by the form $ g_{0 \mu}^{(\mathrm{\rm{GI}, 1)}} = 0 $ and the scalar-vector-tensor decomposition. The obtained result is expressed as

      $X^0 = \frac{1}{a} \int \mathrm{d} \eta \{a \phi^{(1)} \}, \tag{26a} $

      $X^j \!=\! \delta^{ji} \int \mathrm{d} \eta \left\{ \nu_i^{(1)} \!+\! \partial_i b^{(1)} \!+\! \frac{1}{a}\! \int \mathrm{d} \eta' \left\{a \partial_i \phi^{(1)} \right\} \right\} . \tag{26b}$

      Then we rewrite the gauge invariant variables in the form of

      $ \Psi^{(1)} = \psi^{(1)} + \frac{\dot{a}}{a^2} \int \mathrm{d} \eta \{a \phi^{(1)} \}, \tag{27a} $

      $ E^{(1)} = {e}^{(1)} + \int \mathrm{d} \eta \left\{ b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta' \left\{a \phi^{(1)} \right\} \right\}, \tag{27b} $

      $ C_i^{(1)} = c_i^{(1)} - \int \nu_i^{(1)} \mathrm{d} \eta, \tag{27c} $

      $ H_{ij}^{(1)} = h_{ij}^{(1)} . \tag{27d}$

      The above approach has been used to study the scalar cosmological perturbations in Ref. [56]. As suggested in Ref. [48], the gauge invariant synchronous variables are not unique, owing to the indefinite integral in Eq. (26). In this sense, it might be difficult to define observables with the gauge invariant synchronous variables in Eqs. (27a)-(27d).

    • D.   Conversion among different families of gauge invariant first order variables

    • In general, the gauge invariant variables are not limited to the forms in Eqs. (16a)-(16c) and (24a)-(24c). There are other families of gauge invariant variables (see reviews in Refs. [11, 57]). For two different families of the gauge invariant first order variables of metric perturbations $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 1)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} $, the conversion between them can be derived from

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 1)}} - g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} = &(g_{\mu \nu}^{(1)} -{\cal{L}}_{X^A} g_{\mu \nu}^{(0)}) - (g_{\mu \nu}^{(1)} -{\cal{L}}_{X^B} g_{\mu \nu}^{(0)}) \\ =& {\cal{L}}_{(X^B - X^A)} g_{\mu \nu}^{(0)} . \end{aligned} $

      (28)

      If we let $ Z_1^{\mathrm{AB}} \equiv X^A - X^B $, Eq. (28) can be rewritten as

      $ \begin{array}{l} g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 1)}} = g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} -{\cal{L}}_{Z_1^{\mathrm{AB}}} g_{\mu \nu}^{(0)} . \end{array} $

      (29)

      The variable $ Z_1^{\mathrm{AB}} $ was mentioned by Nakamura [54]. It relates two different families of gauge invariant first order variables of metric perturbations. In this study, we further demonstrate that the infinitesimal vector $ Z_1 $ can be expressed as a linear combination of the gauge invariant variables $ {\cal{A}}^{(1)} $, $ {\cal{B}}^{(1)} $ and $ {\cal{C}}^{(1)}_i $:

      $ \begin{array}{l} Z^{\mu, \mathrm{AB}}_1 = \hat{N}_1^{\mu} {\cal{A}}^{(1)} + \hat{N}_2^{\mu} {\cal{B}}^{(1)} + \hat{N}_3^{\mu i} {\cal{C}}_i^{(1)}, \end{array} $

      (30)

      where $ \hat{N}^{\mu}_{1} $, $ \hat{N}^{\mu}_{2} $, and $ \hat{N}^{\mu j}_{3} $ are arbitrary linear operators that are irrelative to any perturbations, and we obtain the gauge invariant variables,

      $ {\cal{A}}^{(1)} \equiv \partial_0 \left( \frac{a^2}{\dot{a}} \psi^{(1)} \right) + a \phi^{(1)}, \tag{31a} $

      $ {\cal{B}}^{(1)} \equiv \partial_0 {e}^{(1)} - b^{(1)} + \frac{a}{\dot{a}} \psi^{(1)}, \tag{31b}$

      $ {\cal{C}}_i^{(1)} \equiv \nu_i^{(1)} - \partial_0 c_i^{(1)} . \tag{31c} $

      The above expressions of $ {\cal{A}}^{(1)} $, $ {\cal{B}}^{(1)} $ and $ {\cal{C}}^{(1)}_i $ can be obtained by making use of Eq. (15).

      The existence of infinite families of gauge invariant variables can also be indicated by the infinite number of choices of $ X^{\mu} $. To be specific, we can extend the expression of Eq. (22) to be

      $ X^0 = \partial_0 {e}^{(1)} - b^{(1)} + Z^0_1, \tag{32a}$

      $ X^j = \delta^{jk} (c_k^{(1)} + \partial_k {e}^{(1)}) + Z^j_1 . \tag{32b} $

      In this way, the gauge invariant metric perturbations could take a general form of

      $ g_{00}^{(\mathrm{\rm{GI}, 1)}} = - 2 a^2 \Phi^{(1)}, \tag{33a}$

      $ g_{0 i}^{(\mathrm{\rm{GI}, 1)}} = a^2 (\partial_i B^{(1)} + V_i^{(1)}), \tag{33b} $

      $ \begin{aligned}[b] g_{ij}^{(\mathrm{\rm{GI}, 1)}} =& a^2 \left(- 2 \Psi^{(1)} \delta_{ij} + 2 \partial_i \partial_j E^{(1)}\right.\\ &\left.+ \partial_i C_j^{(1)} + \partial_j C_i^{(1)} + H_{ij}^{(1)}\right),\end{aligned} \tag{33c} $

      where the gauge invariant variables are defined as

      $ \Phi^{(1)} = \phi^{(1)} - \frac{1}{a} \partial_0 (a (\partial_0 {e}^{(1)} - b^{(1)} + Z^0_1)), \tag{34a}$

      $ B^{(1)} = Z^0_1 - \partial_0 \Delta^{- 1} \partial_j Z^j_1, \tag{34b}$

      $ \Psi^{(1)} = \psi^{(1)} + \frac{\dot{a}}{a} (\partial_0 {e}^{(1)} - b^{(1)} + Z^0), \tag{34c} $

      $ E^{(1)} = \Delta^{- 1} \partial_j Z^j_1, \tag{34d} $

      $ V_i^{(1)} = \nu_i^{(1)} - \partial_0 c_i^{(1)} - \left(\delta_{ik} - \partial_i \Delta^{- 1} \partial_k\right) \partial_0 Z^k_1, \tag{34e}$

      $ C_i^{(1)} = - \left(\delta_{ik} - \partial_i \Delta^{- 1} \partial_k\right) Z^k_1, \tag{34f} $

      $ H_{ij}^{(1)} = h_{ij},\tag{34g} $

      and $ \Delta^{- 1} $ is the inverse Laplacian operator on the background.

      As an example, we consider that $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 1)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} $ are synchronous in Eqs. (27a)-(27d), and Newtonian in Eqs. (23a)-(23d) variables, respectively. In this case, we obtain the explicit expression of $ Z_1^{\mu, \mathrm{AB}} $ as

      $ Z^{0, \mathrm{AB}}_1 = \frac{1}{a} \int {\cal{A}}^{(1)} \mathrm{d} \eta -{\cal{B}}^{(1)}, \tag{35a} $

      $ \begin{aligned}[b] Z^{j, \mathrm{AB}}_1 =& \delta^{jk} \partial_k \int \mathrm{d} \eta \left\{ \frac{1}{a} \int {\cal{A}}^{(1)} \mathrm{d} \eta' \right\} \\ &- \delta^{jk} \partial_k \int {\cal{B}}^{(1)} \mathrm{d} \eta + \delta^{ji} \int {\cal{C}}_i^{(1)} \mathrm{d} \eta . \end{aligned} \tag{35b}$

      Therefore, these two families of gauge invariant variables of metric perturbations are related via the expressions of the linear operators $ \hat{N}_1^{\mu} $, $ \hat{N}_2^{\mu} $, and $ \hat{N}_3^{\mu i} $.

    IV.   GAUGE INVARIANT VARIABLES FOR THE SECOND ORDER COSMOLOGICAL METRIC PERTURBATIONS
    • In this section, the gauge invariant second order variables in cosmology will be derived in the framework of the Lie derivative method and scalar-vector-tensor decomposition. Previous similar studies can be found in Refs. [11, 54]. Furthermore, we will present the conversion formulae among different families of gauge invariant variables. For illustrations, we will consider the gauge invariant metric perturbations that are Newtonian and synchronous, respectively.

    • A.   Gauge transformations of the second order metric perturbations

    • Based on Eq. (3b), the gauge transformation of the second order metric perturbations is presented explicitly as follows

      $ \begin{aligned}[b] \tilde{g}^{(2)}_{00} = &- 2 a^2 \tilde{\phi}^{(2)} \\ =& - 2 a^2 \phi^{(2)} + 2{\cal{L}}_{\xi_1} g_{00}^{(1)} + \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right) g_{00}^{(0)}, \end{aligned}\tag{36a} $

      $ \begin{aligned}[b] \tilde{g}_{0 i}^{(2)} = &a^2 \left(\partial_i \tilde{B}^{(2)} + \tilde{\nu}_i^{(2)}\right) \\ =& g_{00}^{(2)} + 2{\cal{L}}_{\xi_1} g_{0 i}^{(1)} + \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right) g_{0 i}^{(0)} \\ =& a^2 \left(\partial_i \left(b^{(2)} + \Delta^{- 1} \partial^j \Xi_{0 j}\right) + \nu_i^{(2)} + \left(\delta_i^j - \partial_i \Delta^{- 1} \partial^j\right) \Xi_{0 j}\right) \\ \equiv & a^2 \left(\partial_i \left(b^{(2)} + \Delta^{- 1} \partial^j \Xi_{0 j}\right) + \nu_i^{(2)} +{\cal{T}}^j_i \Xi_{0 j}\right), \end{aligned} \tag{36b}$

      $ \begin{aligned}[b] \tilde{g}^{(2)}_{ij} =& a^2 (- 2 \delta_{ij} \tilde{\psi}^{(2)} + 2 \partial_i \partial_j \tilde{e}^{(2)} + \partial_i \tilde{c}^{(2)}_j + \partial_j \tilde{c}_i^{(2)} + \tilde{h}_{ij}^{(2)}) \\ =& g_{ij}^{(2)} + 2{\cal{L}}_{\xi_1} g_{ij}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) g_{ij}^{(0)} \\ =& a^2 \Biggr(- 2 \delta_{ij} \psi^{(2)} + 2 \partial_i \partial_j {e}^{(2)} + \partial_i {c}^{(2)}_j + \partial_j c_i^{(2)} + h_{ij}^{(2)} \\ & + \frac{1}{2} \delta_{ij} \left(\delta^{kl} - \partial^k \Delta^{- 1} \partial^l\right) \Xi_{kl} + \frac{1}{2} \partial_i \partial_j \Delta^{- 1} \left(\left(3 \Delta^{- 1} \partial^k \partial^l\right.\right. \\ &\left.\left.- \delta^{kl}\right) \Xi_{kl}\right) + \partial_j \Delta^{- 1} \partial^l \left(\left(\delta^k_i - \partial_i \Delta^{- 1} \partial^k\right) \Xi_{kl}\right) \\ & + \!\partial_i \Delta^{- 1} \partial^k \left( \left(\delta^l_j \!-\! \partial_j \Delta^{- 1} \partial^l\right) \Xi_{kl}\right)\! +\! \Biggr( \left(\delta^k_i \!-\! \partial_i \Delta^{- 1} \partial^k\right) \left(\delta^l_j \right. \\ &\left. - \partial_j \Delta^{- 1} \partial^l\right) - \frac{1}{2} \left(\delta_{ij} - \partial_i \Delta^{- 1} \partial_j\right) \left(\delta^{kl} - \partial^k \Delta^{- 1} \partial^l\right) \Biggr) \Xi_{kl} \Biggr) \\ =& a^2 \Biggr( - 2 \delta_{ij} \psi^{(2)} + 2 \partial_i \partial_j {e}^{(2)} + \partial_i {c}^{(2)}_j + \partial_j c_i^{(2)} + h_{ij}^{(2)} \\ &+ \frac{1}{2} \delta_{ij} {\cal{T}}^{kl} \Xi_{kl} + \partial_i \partial_j \Delta^{- 1} \left(\left( \partial^k \Delta^{- 1} \partial^l - \frac{1}{2} {\cal{T}}^{kl} \right) \Xi_{kl}\right) \\ &+ \partial_j \Delta^{- 1} \partial^l {\cal{T}}^k_i \Xi_{kl} + \partial_i \Delta^{- 1} \partial^k {\cal{T}}^l_j \Xi_{kl} \\ &+ \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) \Xi_{kl} \Biggr), \end{aligned} \tag{36c}$

      where $ {\cal{T}}^l_j \equiv \delta^l_j - \partial_j \Delta^{- 1} \partial^l $ is a transverse operator, and we define

      $ \Xi_{\mu \nu} \equiv \frac{2{\cal{L}}_{\xi_1} g_{\mu \nu}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) g_{\mu \nu}^{(0)}}{a^2} . \tag{37} $

      In Eqs. (36b) and (36c), we have decomposed the gauge transformation into the scalar, vector, and tensor components. Therefore, the gauge transformation of each component can be rewritten as

      $ \tilde{\phi}^{(2)} = \phi^{(2)} - \frac{1}{2} \Xi_{00}, \tag{38a} $

      $ \tilde{b}^{(2)} = b^{(2)} + \Delta^{- 1} \partial^j \Xi_{0 j}, \tag{38b} $

      $ \tilde{\nu}^{(2)}_i = \nu_i^{(2)} +{\cal{T}}^j_i \Xi_{0 j}, \tag{38c}$

      $ \tilde{\psi}^{(2)} = \psi^{(2)} - \frac{1}{4} {\cal{T}}^{kl} \Xi_{kl}, \tag{38d} $

      $ \tilde{e}^{(2)} = {e}^{(2)} + \frac{1}{2} \Delta^{- 1} \left(\left( \partial^k \Delta^{- 1} \partial^l - \frac{1}{2} {\cal{T}}^{kl} \right) \Xi_{kl}\right), \tag{38e} $

      $ \tilde{c}_j^{(2)} = c_j^{(2)} + \Delta^{- 1} \partial^k {\cal{T}}^l_j \Xi_{kl},\tag{38f} $

      $ \tilde{h}_{ij}^{(2)} = h^{(2)}_{ij} + \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) \Xi_{kl} . \tag{38g} $

      In particular, the second order tensor perturbation is no longer invariant upon the gauge transformation. Based on Eqs. (3b), (4b), and (38), we obtain the gauge invariant second order variables as

      $ \Phi^{(2)} = \phi^{(2)} - \left( \partial_0 + \frac{\dot{a}}{a} \right) Y^0 +{\cal{X}}_{00}, \tag{39a} $

      $ B^{(2)} = b^{(2)} - \partial_0 Y_S + Y^0 - \Delta^{- 1} \partial^j {\cal{X}}_{0 j}, \tag{39b} $

      $ V_i^{(2)} = \nu_i^{(2)} - \delta_{ij} \partial_0 Y^j_T -{\cal{T}}^j_i {\cal{X}}_{0 j},\tag{39c}$

      $ \Psi^{(2)} = \psi^{(2)} + \frac{\dot{a}}{a} Y^0 + \frac{1}{4} {\cal{T}}^{kl} {\cal{X}}_{kl} \tag{39d} $

      $ E^{(2)} \!=\! {e}^{(2)} \!-\! Y_S \!-\! \frac{1}{2} \Delta^{- 1} \left(\left( \partial^k \Delta^{- 1} \partial^l \!-\! \frac{1}{2} {\cal{T}}^{kl} \right) {\cal{X}}_{kl}\right), \tag{39e} $

      $ C_j^{(2)} = c_j^{(2)} - \delta_{jk} Y^k_T - \Delta^{- 1} \partial^k {\cal{T}}^l_j {\cal{X}}_{kl}, \tag{39f} $

      $ H_{ij}^{(2)} = h_{ij}^{(2)} - \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) {\cal{X}}_{kl},\tag{39g} $

      where we define

      $ {\cal{X}}_{\mu \nu} \equiv \frac{2{\cal{L}}_X g_{\mu \nu}^{(1)} -{\cal{L}}_X^2 g_{\mu \nu}^{(0)}}{a^2} , \tag{40}$

      In addition, we decomposed $ Y^i = : Y_T^i + \delta^{ij} \partial_j Y_S $, i.e., into the transverse and longitudinal parts. We infer that $ {\cal{X}}_{\mu\nu} $ depends on $ X^{\mu} $, which determines the gauge invariant first order variables. In Appendix E, we present an explicit expression of $ {\cal{X}}_{\mu \nu} $. As shown in Eq. (39), all the gauge invariant variables are expressed in terms of $ X^{\mu} $ and $ Y^{\mu} $, except that the gauge invariant second order tensor perturbation $ H^{(2)}_{ij} $ solely depends on $ X^{\mu} $. Because the explicit expression of $ X^{\mu} $ has been known, only $ Y^{\mu} $ is undetermined in Eq. (39). Therefore, we will show how to express $ Y^{\mu} $ in terms of the first and second order metric perturbations in the following.

    • B.   Gauge invariant second order Newtonian variables

    • To obtain the gauge invariant Newtonian variables, we derive the expression of $ Y^\mu $ from $ B^{(2)} = E^{(2)} = C_j^{(2)} = 0 $. The explicit expression of $ Y^{\mu} $ is given by

      $ \begin{aligned}[b] Y^0 =& \partial_0 {e}^{(2)} - b^{(2)} + \Delta^{- 1} \partial^j {\cal{X}}_{0 j} - \frac{1}{2} \Biggr( \partial^k \Delta^{- 2} \partial^l \partial_0 \\ &- \frac{1}{2} \Delta^{- 1}{\cal{T}}^{kl}\partial_0 \Biggr) {\cal{X}}_{kl}, \end{aligned}\tag{41a} $

      $ \begin{aligned}[b] Y^i =& \delta^{ij} (\partial_j {e}^{(2)} + c_j^{(2)}) + \Biggr( \frac{1}{4} \partial^i\partial^k \Delta^{- 2} \partial^l \\ &+ \frac{1}{4}\delta^{kl} \partial^i \Delta^{- 1} - \delta^{li} \Delta^{- 1} \partial^k \Biggr) {\cal{X}}_{kl} . \end{aligned}\tag{41b} $

      We determine that $ Y^{\mu} $ depends on a choice of the first order variable $ X^{\mu} $. Therefore, the gauge invariant second order metric perturbations become

      $ g^{(\mathrm{\rm{GI}, 2)}}_{00} = - 2 a^2 \Phi^{(2)}, \tag{42a} $

      $ g_{0 i}^{(\mathrm{\rm{GI}, 2)}} = a^2 V_i^{(2)}, \tag{42b} $

      $ g_{ij}^{(\mathrm{\rm{GI}, 2)}} = - 2 a^2 \delta_{ij} \Psi^{(2)} + a^2 H_{ij}^{(2)}, \tag{42c} $

      where the gauge invariant variables are defined as

      $ \begin{aligned}[b] \Phi^{(2)} = &\phi^{(2)} - \left( \frac{\dot{a}}{a} + \partial_0 \right) \left(\partial_0 {e}^{(2)} - b^{(2)}\right) +{\cal{X}}_{00} \\ &- \left(\frac{\dot{a}}{a} + \partial_0 \right) \left(\Delta^{- 1} \partial^j {\cal{X}}_{0 j} - 3 \Delta^{- 2} \partial^k \partial^l {\cal{X}}_{kl} \right.\\ & \left.+ \delta^{kl} \Delta^{- 1}{\cal{X}}_{kl}\right), \end{aligned} \tag{43a}$

      $ \begin{aligned}[b] \Psi^{(2)} =& \psi^{(2)} + \frac{\dot{a}}{a} (\partial_0 {e}^{(2)} - b^{(2)})\\ &+ \frac{\dot{a}}{a} (\Delta^{- 1} \partial^j {\cal{X}}_{0 j} - 3 \Delta^{- 2} \partial^k \partial^l {\cal{X}}_{kl} \\ & + \delta^{kl} \Delta^{- 1}{\cal{X}}_{kl}) +{\cal{T}}^{kl} {\cal{X}}_{kl}, \end{aligned} \tag{43b}$

      $ V_i^{(2)} = \nu_i^{(2)} - \partial_0 c_i^{(2)} + \Delta^{- 1} \partial^k {\cal{T}}^l_i \partial_0 {\cal{X}}_{kl} -{\cal{T}}^j_i {\cal{X}}_{0 j},\tag{43c} $

      $ H_{ij}^{(2)} = h^{(2)}_{ij} - \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) {\cal{X}}_{kl} . \tag{43d}$

    • C.   Gauge invariant second order synchronous variables

    • For gauge invariant synchronous variables, the $ Y^\mu $ is determined by adopting $ \Phi^{(2)} = 0 $, $ B^{(2)} = 0 $, and $ V^{(2)}_i = 0 $. The explicit expression of the $ Y^\mu $ takes the form of

      $ Y^0 = \frac{1}{a} \int \mathrm{d} \eta \{ a \phi^{(2)} + \frac{1}{2} a{\cal{X}}_{00} \}, \tag{44a} $

      $ \begin{aligned}[b] Y^k =& \delta^{ki} \int \mathrm{d} \eta \left\{ \nu_i^{(2)} + \partial_i b^{(2)} + \frac{1}{a} \int \mathrm{d} \eta \left\{a \partial_i \phi^{(2)} \right\} \right\} \\ &- \int \mathrm{d} \eta \left\{ \delta^{kj} {\cal{X}}_{0 j} - \frac{1}{2 a} \int \mathrm{d} \eta \left\{a \partial^k {\cal{X}}_{00} \right\} \right\} . \end{aligned}\tag{44b} $

      Therefore, we obtain the gauge invariant second order metric perturbations in the form

      $ g_{00}^{(\mathrm{\rm{GI}, 2)}} = 0, \tag{45a} $

      $ g_{0 i}^{(\mathrm{\rm{GI}, 2)}} = 0, \tag{45b} $

      $ \begin{aligned}[b] g_{ij}^{(\mathrm{\rm{GI}, 2)}} = & a^2 \left(- 2 \Psi^{(2)} \delta_{ij} + 2 \partial_i \partial_j E^{(2)} \right.\\ & \left.+ \partial_i C_j^{(2)} + \partial_j C_i^{(2)} + H_{ij}^{(2)}\right),\end{aligned} \tag{45c} $

      where the gauge invariant variables are defined with

      $ \Psi^{(2)} \!=\! \psi^{(2)} \!+\! \frac{\dot{a}}{a^2} \int \mathrm{d} \eta \left\{ a \phi^{(2)} \!+\! \frac{1}{2} a{\cal{X}}_{00} \right\} \!+\!{\cal{T}}^{kl} {\cal{X}}_{kl}, \tag{46a} $

      $ \begin{aligned}[b] E^{(2)} =& {e}^{(2)} - \int \mathrm{d} \eta \left\{ b^{(2)} + \frac{1}{a} \int \mathrm{d} \eta' \left\{a \phi^{(2)} \right\} \right\} \\ & + \int \mathrm{d} \eta \left\{ \Delta^{- 1} \partial^j {\cal{X}}_{0 j} - \frac{1}{2 a} \int \mathrm{d} \eta' \left\{a{\cal{X}}_{00} \right\} \right\} \\ & - \frac{1}{4} (3 \Delta^{- 2} \partial^k \partial^l - \delta^{kl}\Delta^{- 1}) {\cal{X}}_{kl}, \end{aligned} \tag{46b}$

      $ C_j^{(2)}\! \!=\! c_j^{(2)} \!-\!\! \!\int \nu_j^{(2)} \mathrm{d} \eta \!+\!{\cal{T}}^k_j \!\!\!\int\!\! \!{\cal{X}}_{0 k} \mathrm{d} \eta \!-\! \Delta^{- 1} \partial^k {\cal{T}}^l_j {\cal{X}}_{kl}, \tag{46c}$

      $ H_{ij}^{(2)} = h^{(2)}_{ij} - \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) {\cal{X}}_{kl} . \tag{46d} $

      The gauge invariant synchronous variables are different from the gauge invariant Newtonian variables, except the tensor perturbations. In both of the two cases, $ H_{ij}^{(2)} $ is completely determined by the choice of the gauge invariant first order variables.

    • D.   Conversion among different families of gauge invariant second order variables

    • Compared with the first order case in the previous section, a conversion between two different families of gauge invariant second order variables is more complicated, because the gauge invariant second order variables are also dependent on the choice of the gauge invariant first order variables. Specifically, for the two different families of gauge invariant second order metric perturbations $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 2)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 2)}} $, the conversion between them can be derived from

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 2)}} - g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 2)}} =& g_{\mu \nu}^{(2)} - 2{\cal{L}}_{X^A} g_{\mu \nu}^{(1)} - \left({\cal{L}}_{Y^A} -{\cal{L}}_{X^A}^2\right) g_{\mu \nu}^{(0)} \\ &- \left(g_{\mu \nu}^{(2)} - 2{\cal{L}}_{X^B} g_{\mu \nu}^{(1)} - \left({\cal{L}}_{Y^B} -{\cal{L}}_{X^B}^2\right) g_{\mu \nu}^{(0)}\right) \\ = &2{\cal{L}}_{(X^B - X^A)} \left(g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} +{\cal{L}}_{X^B} g_{\mu \nu}^{(0)}\right) \\ &+ \left({\cal{L}}_{(Y^B - Y^A)} -{\cal{L}}_{X^B}^2 -{\cal{L}}_{X^A}^2\right) g_{\mu \nu}^{(0)} \\ =& - 2{\cal{L}}_{(X^A - X^B)} g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} \\ &- \left({\cal{L}}_{(Y^A - Y^B + [X^A, X^B])}\right.\\ &\left.-{\cal{L}}_{(X^A - X^B)}^2\right) g_{\mu \nu}^{(0)} . \end{aligned} $

      (47)

      If we let $ Z_2^{\mathrm{AB}} \equiv Y^A - Y^B + [X^A, X^B] $ and $ Z_1^{\mathrm{AB}} \equiv X^A - X^B $, Eq. (47) can be rewritten as

      $ \begin{array}{l} g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 2)}} \!=\! g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 2)}} \!-\! 2{\cal{L}}_{Z^{\mathrm{AB}}_1} g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} \!-\! ({\cal{L}}_{Z_2^{\mathrm{AB}}} \!-\!{\cal{L}}_{Z^{\mathrm{AB}}_1}^2) g_{\mu \nu}^{(0)} . \end{array} $

      (48)

      This conversion was also mentioned by Nakamura [54] in a different formula. It can be easily verified that both $ Z_1^{\mathrm{AB}} $ and $ Z_2^{\mathrm{AB}} $ are gauge invariant by adopting Eqs. (5a) and (5b). The above formula is generic. It is evident that the gauge invariant second order variables depend on the gauge invariant first order variables, because, in general, we have $ Z_{1}^{\mathrm{AB}}\neq 0 $, i.e., two different families of gauge invariant first order variables. This generic case will be studied later in this section. When we take the same family of the gauge invariant variables at the first order, i.e., $ Z_{1}^{\mathrm{AB}} = 0 $, Eq. (48) can be reduced to a simpler form $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A, 2)}} = g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 2)}}-{\cal{L}}_{Z_2^{\mathrm{AB}}}g_{\mu \nu}^{(0)} $, where $ Z_2^{\mathrm{AB}} = Y^A - Y^B $. For this simple case, the formula is similar to that of Eq. (29).

      We also have infinite families of gauge invariant second order variables. Similar to $ Z_1^{\mathrm{AB}} $ in Eq. (30), the infinitesimal vector $ Z_2^{\mathrm{AB}} $ can also be expressed as a linear combination of the gauge invariant second order variables $ {\cal{A}}^{(2)} $, $ {\cal{B}}^{(2)} $, $ {\cal{C}}^{(2)}_i $ and $ {\cal{D}}_{\sigma \kappa}^{(2)} $,

      $ \begin{aligned}[b]Z_2^{0, \mathrm{AB}} =& \hat{M}_1^0 {\cal{A}}^{(2)} + \hat{M}_2^0 {\cal{B}}^{(2)} + \hat{M}_3^{0 i} {\cal{C}}_i^{(2)} \\& + \hat{M}^{0, \sigma \kappa}_4 {\cal{D}}_{\sigma \kappa}^{(2)}, \end{aligned} \tag{49a}$

      $\begin{aligned}[b] Z_2^{k, \mathrm{AB}} =& \hat{M}_1^k {\cal{A}}^{(2)} + \hat{M}_2^k {\cal{B}}^{(2)} + \hat{M}_3^{ki} {\cal{C}}_i^{(2)} \\&+ \hat{M}^{k, \sigma \kappa}_4 {\cal{D}}_{\sigma \kappa}^{(2)},\end{aligned} \tag{49b} $

      where $ \hat{M}_1^\mu $, $ \hat{M}_2^\mu $, $ \hat{M}_3^{\mu i} $, and $ \hat{M}^{\mu, \sigma \kappa}_4 $ are four arbitrary linear operators that are irrelative to any perturbations, and the gauge-invariant variables are defined as

      $ \begin{aligned}[b] {\cal{A}}^{(2)} \equiv &\partial_0 \left( \frac{a^2}{\dot{a}} \psi^{(2)} \right) + a \phi^{(2)} + \frac{1}{4} \partial_0 \left( \frac{a^2}{\dot{a}} {\cal{T}}^{kl} {\cal{X}}_{kl} \right) \\ &+ \frac{1}{2} a{\cal{X}}_{00}, \end{aligned}\tag{50a} $

      $ \begin{aligned}[b] {\cal{B}}^{(2)} \equiv &\partial_0 {e}^{(2)} - b^{(2)} + \frac{a}{\dot{a}} \psi^{(2)} + \Delta^{- 1} \partial^j {\cal{X}}_{0 j} \\ &+ \! \frac{1}{4} \left( \frac{a}{\dot{a}} {\cal{T}}^{kl} \!-\! 3 \Delta^{- 2} \partial^k \partial^l \partial_0 \!+\! \Delta^{- 1}\delta^{kl}\partial_0 \right) {\cal{X}}_{kl}, \end{aligned} \tag{50b}$

      $ {\cal{C}}_k^{(2)} \equiv \nu_k^{(2)} - \partial_0 c_k^{(2)} + \Delta^{- 1} \partial^i {\cal{T}}^l_k \partial_0 {\cal{X}}_{il} -{\cal{T}}^l_k {\cal{X}}_{0 l}, \tag{50c} $

      $ {\cal{D}}_{\mu \nu}^{(2)} \equiv \frac{1}{a^2}{\left(2{\cal{L}}_{Z_1^{\mathrm{AB}}} g_{\mu \nu}^{(\mathrm{\rm{GI}, B, 1)}} -{\cal{L}}_{Z_1^{\mathrm{AB}}}^2 g_{\mu \nu}^{(0)}\right)} . \tag{50d} $

      Here, we express $ {\cal{A}}^{(2)} $, $ {\cal{B}}^{(2)} $, and $ {\cal{C}}^{(2)}_i $ in terms of $ {\cal{X}}_{\mu \nu} $, which is completely determined by $ X^{B} $. Based on Eqs. (49a) and (49b), we can also obtain a new family of the gauge invariant second order variables by a conversion from a given family of the gauge invariant second order variables. This prediction is similar to that of the first order case.

      We consider three typical cases in the following. In the case that $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A,n)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B,n)}} $ ($ n = 1,2 $) are synchronous and Newtonian, respectively, we infer that $ Z_1^{\mu, \mathrm{AB}} $ takes the form of Eqs. (35a) and (35b), and $ Z_2^{\mu, \mathrm{AB}} $ is expressed as

      $ Z^{0, \mathrm{AB}}_2 = \frac{1}{a} \int {\cal{A}}^{(2)} \mathrm{d} \eta -{\cal{B}}^{(2)} + \frac{1}{2 a} \int \mathrm{d} \eta \left\{a{\cal{D}}_{00}^{(2)} \right\}, \tag{51a}$

      $ \begin{aligned}[b] Z_2^{j, \mathrm{AB}} = &\delta^{jk} \partial_k \int \mathrm{d} \eta \left\{ \frac{1}{a} \int {\cal{A}}^{(2)} \mathrm{d} \eta' \right\} \\ &- \delta^{jk} \partial_k \int {\cal{B}}^{(2)} \mathrm{d} \eta + \delta^{j k} \int {\cal{C}}_k^{(2)} \mathrm{d} \eta \\ &- \!\int \mathrm{d} \eta \left\{ \delta^{j k} {\cal{D}}_{0 k}^{(2)} \!-\! \frac{1}{2 a} \!\int \!\mathrm{d} \eta \left\{a \partial^j {\cal{D}}_{00}^{(2)} \right\} \right\} . \end{aligned} \tag{51b}$

      Since the vectors $ X^{\mu} $ and $ Y^{\mu} $ are independent, we can choose, e.g., the gauge invariant Newtonian variables at the first order, and the gauge invariant synchronous variables at the second order, and vice versa. First, in the case that $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A,1)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B,1)}} $ are Newtonian while $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A,2)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B,2)}} $ are synchronous and Newtonian, respectively, we obtain $ Z^{\mu, \mathrm{AB}}_1 = 0 $ and

      $ Z^{0, \mathrm{AB}}_2 = \frac{1}{a} \int {\cal{A}}^{(2)} \mathrm{d} \eta -{\cal{B}}^{(2)} , \tag{52a}$

      $ \begin{aligned}[b] Z_2^{j, \mathrm{AB}} = &\delta^{jk} \partial_k \int \mathrm{d} \eta \left\{ \frac{1}{a} \int {\cal{A}}^{(2)} \mathrm{d} \eta' \right\} \\ &- \delta^{jk} \partial_k \int {\cal{B}}^{(2)} \mathrm{d} \eta + \delta^{jk} \int {\cal{C}}_k^{(2)} \mathrm{d} \eta . \end{aligned}\tag{52b} $

      Second, in the case that $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A,1)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B,1)}} $ are synchronous and Newtonian, respectively, while both $ g_{\mu \nu}^{(\mathrm{\rm{GI}, A,2)}} $ and $ g_{\mu \nu}^{(\mathrm{\rm{GI}, B,2)}} $ are Newtonian, we deduce that $ Z_1^{\mu, \rm AB} $ takes the same form as Eqs. (35a) and (35b), and $ Z_2^{\mu, \mathrm{AB}} $ turns to be

      $ \begin{aligned}[b] Z^{0,\rm{AB}}_2 = & \Delta^{- 1} \partial^j {\cal{D}}_{0 j}^{(2)} \\ &- \frac{1}{2} \left( \partial^k \Delta^{- 2} \partial^l\partial_0 - \frac{1}{2} \Delta^{- 1}{\cal{T}}^{kl}\partial_0 \right) {\cal{D}}_{kl}^{(2)}, \end{aligned} \tag{53a} $

      $ Z^{i, \rm{AB}}_2\! \!=\! \left(\! \frac{1}{4} \partial^i\partial^k \Delta^{- 2} \partial^l \!\!+\!\! \frac{1}{4}\delta^{kl} \!\Delta^{- 1}\partial^i \!-\! \delta^{li}\! \Delta^{- 1} \partial^k \right)\! {\cal{D}}_{kl}^{(2)} , \tag{53b} $

      which is expressed in terms of the square of the first order metric perturbations only. The above two cases are called the gauge invariant hybrid variables.

    V.   GAUGE INVARIANT EQUATIONS OF MOTION FOR THE SECOND ORDER COSMOLOGICAL PERTURBATIONS
    • In this section, we will derive the equations of motion of the second order cosmological perturbations, which are obtained from the first order scalar perturbations in the gauge invariant framework. For simplicity, we will consider the gauge invariant Newtonian, synchronous, and hybrid variables that have been introduced in previous sections.

    • A.   Gauge invariant energy-momentum tensor up to second order

    • On the matter side, we expand the energy-momentum tensor of the perfect fluid up to the second order, i.e.,

      $ T_{\mu \nu}^{(\mathrm{\rm{GI})}} \equiv T_{\mu \nu}^{(0)} + T_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} + \frac{1}{2} T_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} +{\cal{O}} \left(T_{\mu \nu}^{(\mathrm{\rm{GI}, 3)}}\right), \tag{54} $

      where

      $ T^{(0)}_{\mu \nu} = u_{\mu}^{(0)} u_{\nu}^{(0)} (\rho^{(0)} + P^{(0)}) + g_{\mu \nu}^{(0)} P^{(0)}, \tag{55a}$

      $ \begin{aligned}[b] T^{(\mathrm{\rm{GI}, 1)}}_{\mu \nu} =& u_{\mu}^{(\mathrm{\rm{GI}, 1)}} u_{\nu}^{(0)} \Bigr(\rho^{(0)} + P^{(0)}\Bigr) + u^{(0)}_{\mu} u^{(\mathrm{\rm{GI}, 1)}}_{\nu} \Bigr(\rho^{(0)} + P^{(0)}\Bigr) \\ &+ u_{\mu}^{(0)} u_{\nu}^{(0)} \Bigr(\rho^{(\mathrm{\rm{GI}, 1)}} + P^{(\mathrm{\rm{GI}, 1)}}\Bigr) \\ &+ g_{\mu \nu}^{(0)} P^{(\mathrm{\rm{GI}, 1)}} + g_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} P^{(0)}, \end{aligned}\tag{55b} $

      $ \begin{aligned}[b] T^{(\rm{GI}, 2)}_{\mu \nu} =& u_{\mu}^{(\rm{GI}, 1)} u_{\nu}^{(\rm{GI}, 1)} \Bigr(\rho^{(0)} + P^{(0)}\Bigr) + u_{\mu}^{(\rm{GI}, 1)} u_{\nu}^{(0)} \Bigr(\rho^{(\rm{GI}, 1)} \\ &+ P^{(\rm{GI}, 1)}\Bigr) + u_{\mu}^{(0)} u_{\nu}^{(\rm{GI}, 1)} \Bigr(\rho^{(\rm{GI}, 1)} + P^{(\rm{GI}, 1)}\Bigr)\\ &\times u_{\mu}^{(\mathrm{\rm{GI}, 2)}} u_{\nu}^{(0)} \Bigr(\rho^{(0)} + P^{(0)}\Bigr) + u^{(0)}_{\mu} u^{(\mathrm{\rm{GI}, 2)}}_{\nu} \Bigr(\rho^{(0)} + P^{(0)}\Bigr) \\ &+ u_{\mu}^{(0)} u_{\nu}^{(0)} \Bigr(\rho^{(\mathrm{\rm{GI}, 2)}} + P^{(\mathrm{\rm{GI}, 2)}}\Bigr) + g_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} P^{(\mathrm{\rm{GI}, 1)}}\\ &+ g_{\mu \nu}^{(0)} P^{(\mathrm{\rm{GI}, 2)}} + g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} P^{(0)} . \end{aligned} \tag{55c}$

      Here, $ \rho^{(0)}, P^{(0)} $ denote the background density and pressure, respectively. $ \rho^{(\mathrm{\rm{GI}}, n)} $, $ u_{\mu}^{(\mathrm{\rm{GI}}, n)} $, and $ P^{(\mathrm{\rm{GI}}, n)} $ denote the gauge invariant n-th order density, pressure, and velocity perturbations, respectively. As has been suggested in Eqs. (8a) and (8b), the gauge invariant matter perturbations are formulated as

      $\rho^{(\mathrm{\rm{GI}, 1)}} = \rho^{(1)} -{\cal{L}}_X \rho^{(0)}, \tag{56a} $

      $ P^{(\mathrm{\rm{GI}, 1)}} = P^{(1)} -{\cal{L}}_X P^{(0)}, \tag{56b} $

      $ u_{\mu}^{(\mathrm{\rm{GI}, 1)}} = u_{\mu}^{(1)} -{\cal{L}}_X u_{\mu}^{(0)}, \tag{56c} $

      $ \rho^{(\mathrm{\rm{GI}, 2)}} = \rho^{(2)} - 2{\cal{L}}_X \rho^{(1)} - \left({\cal{L}}_Y -{\cal{L}}_X^2\right) P^{(0)}, \tag{56d} $

      $ P^{(\mathrm{\rm{GI}, 2)}} = P^{(2)} - 2{\cal{L}}_X P^{(1)} - \left({\cal{L}}_Y -{\cal{L}}_X^2\right) P^{(0)}, \tag{56e}$

      $ u_{\mu}^{(\mathrm{\rm{GI}, 2)}} = u^{(2)}_{\mu} - 2{\cal{L}}_X u_{\mu}^{(1)} - \left({\cal{L}}_Y -{\cal{L}}_X^2\right) u^{(0)}_{\mu}. \tag{56f}$

      The velocity field of the perfect fluid can be redefined as $ (\upsilon^i)^{(\mathrm{\rm{GI}}, n)} \equiv a (u^i)^{(\mathrm{\rm{GI}}, n)} $, where $ (u^0)^{(\mathrm{\rm{GI}}, n)} $ could be determined via $ g_{\mu \nu} u^{\mu} u^{\nu} = - 1 $. The equation of state w and the speed of sound $ c_s $ (and $ c_{s}^{(2)} $) are defined as

      $ P^{(0)} = w \rho^{(0)},\tag{57a} $

      $ P^{(\mathrm{\rm{GI}},1)} = c_s^2 \rho^{(\mathrm{\rm{GI}},1)} , \tag{57b}$

      $ P^{(\mathrm{\rm{GI}},2)} = (c_s^{(2)})^2 \rho^{(\mathrm{\rm{GI}},2)}.\tag{57c} $

      All of them are gauge invariant. In principle, one could choose $ c_{s}^{(2)} $ can be taken to be equal to $ c_{s} $ for the adiabatic perturbation with a constant speed of sound [58].

    • B.   Second order cosmological perturbations induced by the first order Newtonian variables

    • As the first step, we study the gauge invariant Newtonian metric perturbations at the first order, including their equations of motion. The gauge invariant metric up to the first order is given by

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI})}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} =& - a^2 \left(1 + 2 \Phi^{(1)}\right) \mathrm{d} \eta^2 + 2 a^2 V_i^{(1)} \mathrm{d} \eta \mathrm{d} x^i \\ &+ a^2 \delta_{ij} \left(1 - 2 \Psi^{(1)}\right) \mathrm{d} x^i \mathrm{d} x^j, \end{aligned} \tag{58}$

      where the gauge invariant first order Newtonian variables have been expressed in Eq. (23a)-(23c). Here, we neglect the first order tensor perturbation, i.e., $ H^{(1)}_{ij} = 0 $. In contrast to the neglect of the scalar or vector perturbations, neglectig $ H^{(1)}_{ij} $ does not violate the gauge invariance that is defined with all the diffeomorphisms (namely, arbitrary $ \xi_1^{\mu} $). In the second order, we will consider that the gauge invariant metric perturbations are Newtonian and synchronous, respectively.

      Based on Eqs. (10a) and (10b), we express the temporal derivative of the conformal Hubble parameter, and the density and velocity perturbations in terms of the gauge invariant first order metric perturbations, i.e.,

      $ \dot{{\cal{H}}} = - \frac{1}{2} (1 + 3 w) {\cal{H}}^2, \tag{59a}$

      $ \rho^{(0)} = \frac{3{\cal{H}}^2}{\kappa a^2}, \tag{59b} $

      $ \rho^{(\mathrm{\rm{GI}, 1)}} = \frac{- 6{\cal{H}}^2 \Phi^{(1)} - 6{\cal{H}} \partial_0 \Psi^{(1)} + 2 \Delta \Psi^{(1)}}{\kappa a^2}, \tag{59c} $

      $ \upsilon_i^{(\mathrm{\rm{GI}, 1)}} \!=\! - V_i^{(1)} \!+\! \frac{\Delta V_i^{(1)} \!-\! 4 \partial_i ({\cal{H}} \Phi^{(1)} \!+\! 4 \partial_0 \Psi^{(1)})}{6 (1 \!+\! w) {\cal{H}}^2}, \tag{59d} $

      where $ \upsilon_i^{(\mathrm{\rm{GI}, 1)}} = \delta_{ij} (\upsilon^i)^{(\mathrm{\rm{GI}, 1)}} $. By substituting the above equations into the spatial part of the first order Einstein field equation, we obtain

      $ \Phi^{(1)} = \Psi^{(1)}, \tag{60a} $

      $ V^{(1)}_i = 0, \tag{60b} $

      and the equation of motion of $ \Phi^{(1)} $, i.e,

      $ \partial_0^2 \Phi^{(1)} + 3 (1 + c_s^2) {\cal{H}} \partial_0 \Phi^{(1)} + 3 (c^2_s - w) {\cal{H}}^2 \Phi^{(1)} - c_s^2 \Delta \Phi^{(1)} = 0. $

      (61)

      Here, we disregard the decaying mode of the first order vector perturbations. In fact, the above results for the first order cosmological perturbations are well known [8].

    • 1.   Gauge invariant second order Newtonian variables
    • At second order, the gauge invariant Newtonian metric perturbations are given by

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} = &2 a^2 \Phi^{(2)} \mathrm{d} \eta^2 + 2 a^2 V_i^{(2)} \mathrm{d} \eta \mathrm{d} x^i \\ &+ a^2 \left(- 2 \delta_{ij} \Psi^{(2)} + H_{ij}^{(2)}\right) \mathrm{d} x^i \mathrm{d} x^j. \end{aligned} $

      (62)

      The explicit expressions of the gauge invariant second order variables $ \Phi^{(2)} $, $ \Psi^{(2)} $, $ V_i^{(2)} $, and $ H_{ij}^{(2)} $ are presented in Eqs. (43a)-(43d) with the expression of $ {\cal{X}}_{\mu \nu}^N $ in Appendix E. Based on the temporal components of the second order Einstein field equations and Eqs. (59a)-(60b), we can express the gauge invariant second order density perturbations in terms of the gauge invariant metric perturbations, i.e.,

      $ \begin{aligned}[b] \rho^{(\rm{GI}, 2)} = &\frac{2}{3 (1 + w) \kappa a^2 {\cal{H}}^2} \Bigr(9 (1 + w) {\cal{H}}^4 \Bigr(4 \Bigr(\Phi^{(1)}\Bigr)^2 - \Phi^{(2)}\Bigr) \\ &- 9 (1 + w) {\cal{H}}^3 \partial_0 \Psi^{(2)} +{\cal{H}}^2 \Bigr(9 (1 + w) \Bigr(\partial_0 \Phi^{(1)}\Bigr)^2 \\ &+ 24 (1 + w) \Phi^{(1)} \Delta \Phi^{(1)} + 3 (1 + w) \Delta \Psi^{(2)}\\ &+ (5 + 9 w) \partial_i \Phi^{(1)} \partial^i \Phi^{(1)}\Bigr) - 8{\cal{H}} \partial_i \partial_0 \Phi^{(1)} \partial^i \Phi^{(1)} \\ &- 4 \partial_i \partial_0\Phi^{(1)} \partial^i \partial_0\Phi^{(1)}\Bigr) . \end{aligned} $

      (63)

      Using Eqs. (59a)-(60b), and substituting the expression of $ \rho^{(\rm{GI}, 2)} $ into the spatial part of the gauge invariant second order Einstein field equations, we can rewrite the second order Einstein field equations in Eq. (10c) to be

      $ \begin{array}{l} {\cal{G}}_{i j} +{\cal{S}}_{i j} = 0 , \end{array} $

      (64)

      where the tensor $ {\cal{G}}_{i j} $ comprises the second order metric perturbations while the tensor $ {\cal{S}}_{i j} $ consists of the square of the first order metric perturbations. Here, their explicit expressions are presented as follows

      $ \begin{aligned}[b] {\cal{S}}_{i j} = &\delta_{i j} \left( \frac{4 (c_s^{(2)})^2}{3 (1 + w)} \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} - 12 \left(c_s^2 + \left(c_s^{(2)}\right)^2 - 2w\right) {\cal{H}}^2 \left(\Phi^{(1)}\right)^2 - 4 \left(5 + 3 c_s^2\right) {\cal{H}} \Phi^{(1)} \partial_0 \Phi^{(1)} \right. \\ &+ \frac{8 (c_s^{(2)})^2}{3 (1 + w) {\cal{H}}} \partial_k \partial_0 \Phi^{(1)} \partial^k \Phi^{(1)} - \left(1 + 3 \left(c_s^{(2)}\right)^2\right) (\partial_0 \Phi^{(1)})^2 - 4 \Phi^{(1)} \partial_0^2 \Phi^{(1)} - 4 \Phi^{(1)} \Delta \Phi^{(1)} - 3 \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} \\ &\left.- 8 \left(c_s^{(2)}\right)^2 \Phi^{(1)}\Delta \Phi^{(1)} - 3 \left(c_s^{(2)}\right)^2 \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} + \frac{4 (c_s^{(2)})^2}{3 (1 + w) {\cal{H}}^2} \partial_k \partial_0 \Phi^{(1)} \partial^k \partial_0 \Phi^{(1)} + 4 c_s^2 \Phi^{(1)} \Delta \Phi^{(1)} \right) \\ &+ 4 \Phi^{(1)} \partial_i \partial_j \Phi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}} \left(\partial_i \partial_0 \Phi^{(1)} \partial_j \Phi^{(1)} + \partial_i \Phi^{(1)} \partial_0 \partial_j \Phi^{(1)}\right) \\ &+ \left( 2 - \frac{4}{3 (1 + w)} \right) \partial_i \Phi^{(1)} \partial_j \Phi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Phi^{(1)} \partial_i \partial_0 \Phi^{(1)}, \end{aligned} $

      (65)

      $ \begin{aligned}[b] {\cal{G}}_{i j} = & \frac{1}{4} \partial_0^2 H_{i j}^{(2)} + \frac{1}{2} {\cal{H}} \partial_0 H_{i j}^{(2)} - \frac{1}{4} \Delta H_{i j}^{(2)} + \delta_{i j} \left( 3 \left(\left(c_s^{(2)}\right)^2 - w\right) {\cal{H}}^2 \Phi^{(2)} +{\cal{H}} \partial_0 \Phi^{(2)} + \frac{1}{2} \Delta \Phi^{(2)} \right.\\ & + \partial_0^2 \Psi^{(2)} + \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)}\left.- \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} \right) - \frac{1}{2} {\cal{H}} \partial_i V_j^{(2)} \\ & - \frac{1}{4} \partial_i \partial_0 V_j^{(2)} - \frac{1}{2} {\cal{H}} \partial_j V_i^{(2)}- \frac{1}{4} \partial_j \partial_0 V_i^{(2)} - \frac{1}{2} \partial_i \partial_j \Phi^{(2)} + \frac{1}{2} \partial_i \partial_j \Psi^{(2)} . \end{aligned} $

      (66)

      We can decompose Eq. (64) into the tensor, vector, and scalar components. For illustrations, we decompose $ {\cal{G}}_{i j} $ as a first step. The decomposition of $ {\cal{G}}_{i j} $ is explicitly given by

      $ \Lambda_{k l}^{i j} {\cal{G}}_{i j} = \frac{1}{4} \partial_0^2 H_{k l}^{(2)} + \frac{1}{2} {\cal{H}} \partial_0 H_{k l}^{(2)} - \frac{1}{4} \Delta H_{k l}^{(2)} ,\tag{67a} $

      $ \Delta^{- 1} \partial^l {\cal{T}}^k_i {\cal{G}}_{k l} = \frac{1}{4} (\partial_0 + 2{\cal{H}}) V_i^{(2)} , \tag{67b} $

      $ \frac{1}{2} \Delta^{- 1} \left( \partial^k \Delta^{- 1} \partial^l - \frac{1}{2} {\cal{T}}^{k l} \right) {\cal{G}}_{k l} = \frac{1}{4} \left(\Psi^{(2)} - \Phi^{(2)}\right), \tag{67c} $

      $ \begin{aligned}[b] \frac{1}{4} {\cal{T}}^{i j} {\cal{G}}_{i j} = & \frac{1}{2} \left( 3 \left(\left(c_s^{(2)}\right)^2 - w\right) {\cal{H}}^2 \Phi^{(2)} +{\cal{H}} \partial_0 \Phi^{(2)} + \frac{1}{2} \Delta \Phi^{(2)}\right. \\ &+ \partial_0^2 \Psi^{(2)} + \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} \\ &\left.- \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} \right) , \end{aligned} \tag{67d}$

      where $ \Lambda^{kl}_{ij} \equiv {\cal{T}}^k_i {\cal{T}}^l_j - ({1}/{2}) {\cal{T}}_{ij} {\cal{T}}^{kl} $. The equations of motion of the gauge invariant second order cosmological perturbations can be written as

      $\partial_0^2 H_{ij}^{(2)} + 2{\cal{H}} \partial_0 H_{ij}^{(2)} - \Delta H_{ij}^{(2)} = - 4 \Lambda^{kl}_{ij} {\cal{S}}_{kl}, \tag{68a}$

      $ (\partial_0 + 2{\cal{H}}) V_i^{(2)} = - 4 \Delta^{- 1} {\cal{T}}^k_i \partial^l {\cal{S}}_{k l} , \tag{68b}$

      $ \Psi^{(2)} - \Phi^{(2)} = - 2 \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j} , \tag{68c} $

      $ \begin{aligned}[b] \partial_0^2 \Psi^{(2)} +& \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} - \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} \\ +& 3 \left(\left(c_s^{(2)}\right)^2 - w\right) {\cal{H}}^2 \Phi^{(2)} +{\cal{H}} \partial_0 \Phi^{(2)} \\ +& \frac{1}{2} \Delta \Phi^{(2)}= - \frac{1}{2} {\cal{T}}^{i j} {\cal{S}}_{i j} , \end{aligned} \tag{68d}$

      where we have the following expressions of $ {\cal{S}}_{i j} $,

      $ \begin{aligned}[b] \Lambda^{i j}_{k l} {\cal{S}}_{i j} =& \Lambda^{i j}_{k l} \left( \left( 2 + \frac{4}{3 (1 + w)} \right) \Phi^{(1)} \partial_i \partial_j \Phi^{(1)} \right.\\ &+ \frac{8}{3 (1 + w) {\cal{H}}} \Phi^{(1)} \partial_0 \partial_i \partial_j \Phi^{(1)} \\ &\left.- \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Phi^{(1)} \partial_j \partial_0 \Phi^{(1)} \right) , \end{aligned} $

      (69)

      $ \begin{aligned}[b] \Delta^{- 1} {\cal{T}}^k_i \partial^l {\cal{S}}_{k l} =& \Delta^{- 1} {\cal{T}}^k_i \partial^l \left( \left( 2 - \frac{4}{3 (1 + w)} \!\right) \partial_k \Phi^{(1)} \partial_l \Phi^{(1)} \right.- \frac{4}{3 (1 + w) {\cal{H}}} (\partial_k \partial_0 \Phi^{(1)} \partial_l \Phi^{(1)} + \partial_k \Phi^{(1)} \partial_0 \partial_l \Phi^{(1)}) \\ &+ 4 \Phi^{(1)} \partial_k \partial_l \Phi^{(1)} \left.- \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_k \partial_0 \Phi^{(1)} \partial_l \partial_0 \Phi^{(1)} \right) , \end{aligned} $

      (70)

      $ \begin{aligned}[b] \frac{1}{4} {\cal{T}}^{i j} {\cal{S}}_{i j} =& \frac{1}{2} \left( \frac{4 \left(c_s^{(2)}\right)^2}{3 (1 \!+\! w)} \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} \!-\! 12 \left(c_s^2\!+\! \left(c_s^{(2)}\right)^2 \!-\! 2w\right) {\cal{H}}^2 \left(\Phi^{(1)}\right)^2 \!-\! 4 \left(5 \!+\! 3 c_s^2\right) {\cal{H}} \Phi^{(1)} \partial_0 \Phi^{(1)} \!-\! 8 \left(c_s^{(2)}\right)^2 \Phi^{(1)}\Delta \Phi^{(1)}\right. \\ &+ \frac{8 (c_s^{(2)})^2}{3 (1 + w) {\cal{H}}} \partial_k \partial_0 \Phi^{(1)} \partial^k \Phi^{(1)} - \left(1 + 3 \left(c_s^{(2)}\right)^2\right) \left(\partial_0 \Phi^{(1)}\right)^2 - 4 \Phi^{(1)} \partial_0 ^2\Phi^{(1)} - 4 \Phi^{(1)} \Delta \Phi^{(1)} - 3 \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} \\ &\left.- 3 \left(c_s^{(2)}\right)^2 \partial_k \Phi^{(1)} \partial^k \Phi^{(1)} + \frac{4 \left(c_s^{(2)}\right)^2}{3 (1 + w) {\cal{H}}^2} \partial_k \partial_0 \Phi^{(1)} \partial^k \partial_0 \Phi^{(1)} + 4 c_s^2 \Phi^{(1)} \Delta \Phi^{(1)} \right)\\ & + \frac{1}{4} {\cal{T}}^{i j} \left( \left( 2 + \frac{4}{3 (1 + w)} \right) \Phi^{(1)} \partial_i \partial_j \Phi^{(1)} + \frac{8}{3 (1 + w) {\cal{H}}} \Phi^{(1)} \partial_0 \partial_i \partial_j \Phi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Phi^{(1)} \partial_j \partial_0 \Phi^{(1)} \right) , \end{aligned} $

      (71)

      $ \begin{aligned}[b] \frac{1}{2} \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j} =& \frac{1}{2} \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) \left( \left( 2 - \frac{4}{3 (1 + w)} \right) \partial_i \Phi^{(1)} \partial_j \Phi^{(1)} + 4 \Phi^{(1)} \partial_i \partial_j \Phi^{(1)} \right. \\ &- \frac{4}{3 (1 + w) {\cal{H}}} \left(\partial_i \partial_0 \Phi^{(1)} \partial_j \Phi^{(1)} + \partial_i \Phi^{(1)} \partial_0 \partial_j \Phi^{(1)}\right) \left. - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Phi^{(1)} \partial_j \partial_0 \Phi^{(1)} \right) . \end{aligned} $

      (72)

      The above equations of motion can be derived in an easier approach, by adopting the scalar-vector-tensor decomposition operators onto both sides of Eq. (64). These operators are defined on the background, and are thus gauge invariant. Therefore, the equations of motion in Eqs. (68a)-(68d) should also be gauge invariant. Here, the operator $ \Lambda^{kl}_{ij} $ is defined in configuration space. In Appendix F, we express $ \Lambda^{kl}_{ij} $ in momentum space and demonstrate that it can be rewritten in terms of the polarization tensors.

    • 2.   Gauge invariant second order synchronous variables
    • At second order, the gauge invariant synchronous metric perturbations are given by

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} =& a^2 \left(- 2 \delta_{ij} \Psi^{(2)} + 2 \partial_i \partial_j E^{(2)}\right. \\ &\left.+ \partial_i C_j^{(2)} + \partial_j C_i^{(2)} + H_{ij}^{(2)}\right) \mathrm{d} x^i \mathrm{d} x^j . \end{aligned} $

      (73)

      The explicit expressions of the gauge invariant second order variables, $ E^{(2)} $, $ \Psi^{(2)} $, $ C_i^{(2)} $, and $ H_{ij}^{(2)} $ are presented in Eqs. (46a)-(46d) with the expression of $ {\cal{X}}_{\mu \nu}^N $ in Appendix E. Based on the temporal components of the second order Einstein field equations and Eqs. (59a)-(60b), the gauge invariant second order density perturbations in terms of the metric perturbations is given by

      $ \begin{aligned}[b] \rho^{(\rm{GI}, 2)} = &\frac{2}{3 (1 + w) \kappa a^2 {\cal{H}}^2} (36 (1 + w) {\cal{H}}^4 (\Phi^{(1)})^2 \\ &- 3 (1 + w) {\cal{H}}^3 (3 \partial_0 \Psi^{(2)} - \Delta \partial_0 E^{(2)}) \\ &+{\cal{H}}^2 (9 (1 + w) (\partial_0 \Phi^{(1)})^2 + 24 (1 + w) \Phi^{(1)} \Delta \Phi^{(1)} \\ &+ 3 (1 + w) \Delta \Psi^{(2)} + (5 + 9 w) \partial_i \Phi^{(1)} \partial^i \Phi^{(1)}) \\ &- 8{\cal{H}} \partial_i \partial_0 \Phi^{(1)} \partial^i \Phi^{(1)} - 4 \partial_i \partial_0\Phi^{(1)} \partial^i \partial_0\Phi^{(1)}) . \end{aligned} $

      (74)

      Using Eqs. (59a)-(60b) and substituting the expression of $ \rho^{(\rm{GI}, 2)} $ into the spatial parts of the gauge invariant second order Einstein field equations, we can also rewrite the second order Einstein field equations in the terms of the tensors $ {\cal{G}}_{i j} $ and $ {\cal{S}}_{i j} $, i.e.,

      $ \begin{array}{l} {\cal{G}}_{i j} +{\cal{S}}_{i j} = 0, \end{array} $

      (75)

      where the $ {\cal{S}}_{i j} $ takes the same form expressed in Eq. (65), and

      $ \begin{aligned}[b] {\cal{G}}_{i j} =& \frac{1}{4} \partial_0^2 H_{i j}^{(2)} + \frac{1}{2} {\cal{H}} \partial_0 H_{i j}^{(2)} - \frac{1}{4} \Delta H_{i j}^{(2)} + \delta_{i j} \left( \partial_0^2 \Psi^{(2)} + \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} - \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} - \Delta \left( \frac{1}{2} \left(\partial_0 + 2{\cal{H}}\right) \partial_{\eta} E^{(2)}\right.\right.\\ & + \left(c_s^{(2)}\right)^2 {\cal{H}} \partial_{\eta} E^{(2)} \Biggr) \Biggr) + \frac{1}{2} {\cal{H}} \partial_i \partial_0 C_j^{(2)} + \frac{1}{4} \partial_i \partial_0^2 C_j^{(2)} + \frac{1}{2} {\cal{H}} \partial_j \partial_0 C_i^{(2)}+ \frac{1}{4} \partial_j \partial_0^2 C_i^{(2)} +{\cal{H}} \partial_i \partial_j \partial_0 E^{(2)} + \frac{1}{2} \partial_j \partial_j \partial_0^2 E^{(2)} + \frac{1}{2} \partial_i \partial_j \Psi^{(2)} . \end{aligned} $

      (76)

      By adopting the scalar-vector-tensor decomposition to Eq. (75), the equations of motion of the gauge invariant second order cosmological perturbations are obtained as

      $\partial_0^2 H_{i j}^{(2)} + 2{\cal{H}} \partial_0 H_{ i j}^{(2)} - \Delta H_{i j}^{(2)} = - 4 \Lambda^{kl}_{ij} {\cal{S}}_{kl}, \tag{77a} $

      $ (\partial_0 + 2{\cal{H}}) \partial_0 C_i^{(2)} = - 4 \Delta^{- 1} {\cal{T}}^k_i \partial^l {\cal{S}}_{k l}, \tag{77b} $

      $ (\partial_0 + 2{\cal{H}}) \partial_0 E^{(2)} + \Psi^{(2)} = - 2 \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j}, \tag{77c}$

      $ \begin{aligned}[b] \partial_0^2 \Psi^{(2)} &+ \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} - \frac{1}{2} (1 + 2 (c_s^{(2)})^2) \Delta \Psi^{(2)} \\ &- \Delta \left( \frac{1}{2} (\partial_0 + 2{\cal{H}}) \partial_0 E^{(2)} + \left(c_s^{(2)}\right)^2 {\cal{H}} \partial_0 E^{(2)} \right) \\ =& - \frac{1}{2} {\cal{T}}^{i j} {\cal{S}}_{i j} . \end{aligned}\tag{77d} $

      The right hand sides of the above equations are the same as those in Eqs. (69)-(72), because the source term $ {\cal{S}}_{i j} $ is uniquely determined by the gauge invariant first order Newtonian variables. In particular, the equation of motion of the tensor perturbation $ H_{\mu \nu}^{(\rm{GI},2)} $ in Eq. (77a) is as same as that in Eq. (68a). For the vector perturbation, the evolution of $ \partial_0 C_i^{(2)} $ in Eq. (77b) is demonstrated to be the same as that of $ V_i^{(2)} $ in Eq. (68b).

    • C.   Second order cosmological perturbations induced by the first order synchronous variables

    • In the subsection, we turn to discuss the gauge invariant synchronous metric perturbations in the first order, and their equations of motion. For the second order equations of motion, we will also consider the gauge invariant variables that are Newtonian and synchronous, respectively.

      The gauge invariant synchronous metric up to the first order is given by

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI})}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} = &- a^2 \mathrm{d} \eta^2 + a^2 (\delta_{ij} (1 - 2 \Psi^{(1)}) \\ &+ 2 \partial_i \partial_j E^{(1)} + \partial_i C_j^{(1)} + \partial_j C_i^{(1)}) \mathrm{d} x^i \mathrm{d} x^j , \end{aligned} \tag{78}$

      where we also neglect the first order tensor perturbation, i.e., $ H_{i j}^{(1)} = 0 $. Here, the gauge invariant first order variables $ \Psi^{(1)} $, $ E^{(1)} $, and $ C^{(1)}_i $ have been expressed in Eqs. (27a)-(27d). Based on the first-order Einstein field equations and Eqs. (59a) and (59b), the density and velocity perturbations are expressed in terms of the gauge invariant first order metric perturbations, i.e.,

      $ \rho^{(\mathrm{\rm{GI}, 1)}} = \frac{- 6{\cal{H}} \partial_0 \Psi^{(1)} + 2{\cal{H}} \Delta \partial_0 E^{(1)} + 2 \Delta \Psi^{(1)}}{\kappa a^2}, \tag{79a}$

      $ \upsilon_i^{(\mathrm{\rm{GI}, 1)}} = - \frac{\Delta \partial_0 C_i^{(1)} + 4 \partial_i \partial_0 \Psi^{(1)}}{6 (1 + w) {\cal{H}}^2} . \tag{79b} $

      By substituting Eqs. (59a) and (59b), including the equations above, into the spatial parts of the first-order Einstein field equations, we can obtain

      $ C_i^{(1)} = 0, \tag{80}$

      and the equations of motion for $ E^{(1)} $ and $ \Psi^{(1)} $, i.e.,

      $ \Psi^{(1)} + 2{\cal{H}} \partial_0 E^{(1)} + \partial_0^2 E^{(1)} = 0, \tag{81a} $

      $ \partial_0^2 \Psi^{(1)} +{\cal{H}} (2 + 3 c_s^2) \partial_0 \Psi^{(1)} - c_s^2 \Delta \Psi^{(1)} -{\cal{H}}c_s^2 \Delta \partial_0 E^{(1)} = 0. \tag{81b}$

      Here, we also disregard the decaying mode of the first order vector perturbations.

    • 1.   Gauge invariant second order Newtonian variables
    • In the second order, the gauge invariant Newtonian metric perturbations are given by

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} =& 2 a^2 \Phi^{(2)} \mathrm{d} \eta^2 + 2 a^2 V_i^{(2)} \mathrm{d} \eta \mathrm{d} x^i \\ & + a^2 (- 2 \delta_{ij} \Psi^{(2)} + H_{ij}^{(2)}) \mathrm{d} x^i \mathrm{d} x^j, \end{aligned} $

      (82)

      The explicit expressions of the gauge invariant second order variables, $ \Phi^{(2)} $, $ \Psi^{(2)} $, $ V_i^{(2)} $ and $ H_{ij}^{(2)} $ are presented in Eqs. (43a)-(43d) with the expression of $ {\cal{X}}_{\mu \nu}^S $ in Appendix E. Based on the temporal components of the second order Einstein field equations and Eqs. (59a), (59b), (79a)-(80), we can express the gauge invariant second order density perturbation in terms of metric perturbations, namely,

      $ \begin{aligned}[b] \rho^{(\rm{GI}, 2)} =& - \frac{1}{3 (1 + w) \kappa a^2 {\cal{H}}^2} \Big(18 (1 + w) {\cal{H}}^4 \Phi^{(2)} + 8 \partial_i \partial_0 \Psi^{(1)} \partial^i \partial_0 \Psi^{(1)} + 6 (1 + w) {\cal{H}}^3 (3 \partial_0 \Psi^{(2)} - 4 \partial_0 \Psi^{(1)} \Delta E^{(1)} \\ &+ 4 \Psi^{(1)} (3 \partial_0 \Psi^{(1)} - \Delta \partial_0 {e}^{(1)}) + 4 \partial_i \partial_j \partial_0 E^{(1)} \partial^i \partial^j E^{(1)}) - 3 (1 + w) {\cal{H}}^2 (6 (\partial_0 \Psi^{(1)})^2-4 \partial_0 \Psi^{(1)} \Delta \partial_0 E^{(1)} \\ & + 16 \Psi^{(1)} \Delta \Psi^{(1)} + 2 \Delta \Psi^{(2)} + 6 \partial_i \Psi^{(1)} \partial^i \Psi^{(1)} - 4 \partial_i \Delta E^{(1)} \partial^i \Psi^{(1)} + \Delta \partial_0 E^{(1)} \Delta \partial_0 E^{(1)} - 4 \Delta E^{(1)} \Delta \Psi^{(1)} \\ &- \partial_i \Delta E^{(1)} \partial^i \Delta E^{(1)} - 4 \partial_i \partial_j \Psi^{(1)} \partial^i \partial^j E^{(1)} - \partial_i \partial_j \partial_0 E^{(1)} \partial^i \partial^j \partial_0 E^{(1)} + \partial_k \partial_i \partial_j E^{(1)} \partial^i \partial^j \partial^k E^{(1)})\Big). \end{aligned} $

      (83)

      Using Eqs. (59a), (59b), and (79a)-(80), and substituting the expression of $ \rho^{(\rm{GI}, 2)} $ into the spatial parts of the gauge invariant second order Einstein field equations, we can also obtain

      $ \begin{array}{l} {\cal{G}}_{i j} +{\cal{S}}_{i j} = 0, \end{array} $

      (84)

      where the form of $ {\cal{G}}_{i j} $ is the same as in Eq. (66), and

      $ \begin{aligned}[b] {\cal{S}}_{i j} = &\delta_{i j} \Biggr( - 12 \left(c_s^2 - \left(c_s^{(2)}\right)^2\right) {\cal{H}} \Psi^{(1)} \partial_0 \Psi^{(1)} + \left( 1 - 3 \left( c_s^{(2)} \right)^2 \right) \left(\partial_0 \Psi^{(1)}\right)^2 - 2 \left(1 - 2 c_s^2 + 4 \left(c_s^{(2)}\right)^2\right) \Psi^{(1)} \Delta \Psi^{(1)} \\ &- \left(2 + 3 \left(c_s^{(2)}\right)^2\right) \partial_k \Psi^{(1)} \partial^k \Psi^{(1)} + \frac{4 \left(c_s^{(2)}\right)^2}{3 (1 + w) {\cal{H}}^2} \partial_k \partial_0 \Psi^{(1)} \partial^k \partial_0 \Psi^{(1)} - \left(1 + \left(c_s^{(2)}\right)^2\right) \Delta \partial_0 E^{(1)} \Delta \partial_0 E^{(1)}\\ &+ \left(1 + \left(c_s^{(2)}\right)^2\right) \partial_k \Delta E^{(1)} \partial^k \Delta E^{(1)} + 8 \left(1 + \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_k \partial_l \partial_0 E^{(1)} \partial^k \partial^l E^{(1)} + 4 \partial_k \partial_l \partial_0^2 E^{(1)} \partial^k \partial^l E^{(1)}\\ &+ \left(3 + \left(c_s^{(2)}\right)^2\right) \partial_k \partial_l \partial_0 E^{(1)} \partial^k \partial^l \partial_0 E^{(1)} - \left(1 + \left(c_s^{(2)}\right)^2\right) \partial_m \partial_l \partial_k E^{(1)} \partial^m \partial^l \partial^k E^{(1)} \\ &- 2 \partial_0^2 \Psi^{(1)} \Delta E^{(1)} - \left(1 - 2 \left(c_s^{(2)}\right)^2\right) \partial_0 \Psi^{(1)} \Delta \partial_0 E^{(1)} + \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \partial_k \Delta E^{(1)} \partial^k \Psi^{(1)} \\ &+ 2 \left(1 + \left(c_s^{(2)}\right)^2\right) \Delta E^{(1)} \Delta \Psi^{(1)} + 2 \left(c_s^{(2)}\right)^2 \partial_k \partial_l \Psi^{(1)} \partial^k \partial^l E^{(1)} \\ &- 2 - 4{\cal{H}} \Biggr( \partial_0 \Psi^{(1)} \left( 1 + \left(c_s^{(2)}\right)^2 \right) \Delta E^{(1)} + \Psi^{(1)} \left(\left(c_s^{(2)}\right)^2 - c_s^2\right) \Delta \partial_0 E^{(1)} \Biggr) \Biggr)\\ &+ 3 \partial_i \Psi^{(1)} \partial_j \Psi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Psi^{(1)} \partial_j \partial_0 \Psi^{(1)} + 2 \Psi^{(1)} \partial_i \partial_j \Psi^{(1)}\\ &- 2 \partial_k \partial_j \partial_0 E^{(1)} \partial^k \partial_i \partial_0 E^{(1)} - \partial_k \Delta E^{(1)} \partial^k \partial_i \partial_j E^{(1)} + \partial_k \partial_l \partial_j E^{(1)} \partial^k \partial^l \partial_i E^{(1)}\\ &- 2 \Delta \partial_0^2 E^{(1)} \partial_i \partial_j E^{(1)} - 4{\cal{H}} \left(1 + c_s^2\right) \Delta \partial_0 E^{(1)} \partial_i \partial_j E^{(1)} + \Delta \partial_0 E^{(1)} \partial_i \partial_j \partial_0 E^{(1)}\\ &- \partial_k \partial_i \partial_j E^{(1)} \partial^k \Psi^{(1)} + 2 \partial_k \partial_j \Psi^{(1)} \partial^k \partial_i E^{(1)} + 2 \partial_k \partial_i \Psi^{(1)} \partial^k \partial_j E^{(1)} + 12 (1 + c_s^2) {\cal{H}} \partial_0 \Psi^{(1)} \partial_i \partial_j E^{(1)}\nonumber\\ &+ 6 \partial_0^2 \Psi^{(1)} \partial_i \partial_j E^{(1)} - 4 \left(1 + c_s^2\right) \Delta \Psi^{(1)} \partial_i \partial_j E^{(1)} + \partial_0 \Psi^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - 2 \Delta E^{(1)} \partial_i \partial_j \Psi^{(1)} . \end{aligned} \tag{85}$

      Following the approach in the previous subsection, the equations of motion of the gauge invariant second order cosmological perturbations can be written as

      $ \partial_0^2 H_{ij}^{(2)} + 2{\cal{H}} \partial_0 H_{ij}^{(2)} - \Delta H_{ij}^{(2)} = - 4 \Lambda^{kl}_{ij} {\cal{S}}_{kl}, \tag{86a}$

      $ (\partial_0 + 2{\cal{H}}) V_i^{(2)} = - 4 \Delta^{- 1} {\cal{T}}^k_i \partial^l {\cal{S}}_{k l} , \tag{86b}$

      $ \Psi^{(2)} - \Phi^{(2)} = - 2 \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j} , \tag{86c}$

      $ \begin{aligned}[b] \partial_0^2 \Psi^{(2)}& + \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} \\&- \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} \\ &+ 3 \left(\left(c_s^{(2)}\right)^2 - w\right) {\cal{H}}^2 \Phi^{(2)} \\&+{\cal{H}} \partial_0 \Phi^{(2)} + \frac{1}{2} \Delta \Phi^{(2)} \\ =& - \frac{1}{2} {\cal{T}}^{i j} {\cal{S}}_{i j} , \end{aligned}\tag{86d} $

      where

      $ \begin{aligned}[b] \Lambda^{i j}_{k l} {\cal{S}}_{i j} = &\Lambda^{i j}_{k l} \Biggr( -\Psi^{(1)} \partial_i \partial_j \Psi^{(1)} + \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_0 \Psi^{(1)} \partial_i \partial_j \partial_0 \Psi^{(1)}+ \Bigr(- 5 \partial^m \Psi^{(1)} \partial_m - 6 c_s^2 {\cal{H}} \partial_0 \Psi^{(1)} \\ &- 2 \left(1 - c_s^2\right) \Delta \Psi^{(1)} + \partial_0 \Psi^{(1)} \partial_0- 2 \Psi^{(1)} \Delta\Bigr) \partial_i \partial_j E^{(1)} + 2 \partial_m \partial_0 E^{(1)} \partial^m \partial_i \partial_j \partial_0 E^{(1)}+ \Delta \partial_0 E^{(1)} \partial_i \partial_j \partial_0 E^{(1)} \\ &- \partial_m \Delta E^{(1)} \partial^m \partial_i \partial_j E^{(1)}- \partial_n \partial_m E^{(1)} \partial^n \partial^m \partial_i \partial_j E^{(1)} + 2 c_s^2 {\cal{H}} \Delta \partial_0 E^{(1)} \partial_i \partial_j E^{(1)} \Biggr) , \end{aligned} \qquad \qquad \quad $

      (87)

      $ \begin{aligned}[b] \Delta^{- 1} {\cal{T}}^i_m \partial^j {\cal{S}}_{i j} = & \Delta^{- 1} {\cal{T}}^i_m \partial^j \Biggr( 3 \partial_i \Psi^{(1)} \partial_j \Psi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Psi^{(1)} \partial_j \partial_0 \Psi^{(1)} + 2 \Psi^{(1)} \partial_i \partial_j \Psi^{(1)} - 2 \partial_k \partial_j \partial_0 E^{(1)} \partial^k \partial_i \partial_0 E^{(1)} \\ &- \partial_k \Delta E^{(1)} \partial^k \partial_i \partial_j E^{(1)} + \partial_k \partial_l \partial_j E^{(1)} \partial^k \partial^l \partial_i E^{(1)}- 2 \Delta \partial_0^2 E^{(1)} \partial_i \partial_j E^{(1)} - 4{\cal{H}} \left(1 + c_s^2\right) \Delta \partial_0 E^{(1)} \partial_i \partial_j E^{(1)} \\ & + \Delta \partial_0 E^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - \partial_k \partial_i \partial_j E^{(1)} \partial^k \Psi^{(1)}+ 2 \partial_k \partial_j \Psi^{(1)} \partial^k \partial_i E^{(1)} + 2 \partial_k \partial_i \Psi^{(1)} \partial^k \partial_j E^{(1)} + 12 (1 \\ &+ c_s^2) {\cal{H}} \partial_0 \Psi^{(1)} \partial_i \partial_j E^{(1)} + 6 \partial_0^2 \Psi^{(1)} \partial_i \partial_j E^{(1)} - 4 \left(1 + c_s^2\right) \Delta \Psi^{(1)} \partial_i \partial_j E^{(1)} + \partial_0 \Psi^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - 2 \Delta E^{(1)} \partial_i \partial_j \Psi^{(1)} \Biggr) , \end{aligned} $

      (88)

      $ \begin{aligned}[b] \frac{1}{4} {\cal{T}}^{i j} {\cal{S}}_{i j} = & \frac{1}{2} \Biggr( - 12 \left(c_s^2 - \left(c_s^{(2)}\right)^2\right) {\cal{H}} \Psi^{(1)} \partial_0 \Psi^{(1)} + \left( 1 - 3 \left( c_s^{(2)} \right)^2 \right) \left(\partial_0 \Psi^{(1)}\right)^2 - 2 \left(1 - 2 c_s^2 + 4 \left(c_s^{(2)}\right)^2\right) \Psi^{(1)} \Delta \Psi^{(1)} \\ &- \left(2 + 3 \left(c_s^{(2)}\right)^2\right) \partial_k \Psi^{(1)} \partial^k \Psi^{(1)} + \frac{4 (c_s^{(2)})^2}{3 (1 + w) {\cal{H}}^2} \partial_k \partial_0 \Psi^{(1)} \partial^k \partial_0 \Psi^{(1)} - \left(1 + \left(c_s^{(2)}\right)^2\right) \Delta \partial_0 E^{(1)} \Delta \partial_0 E^{(1)} \\ &+ \left(1 + \left(c_s^{(2)}\right)^2\right) \partial_k \Delta E^{(1)} \partial^k \Delta E^{(1)} + 8 \left(1 + \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_k \partial_l \partial_0 E^{(1)} \partial^k \partial^l E^{(1)} + 4 \partial_k \partial_l \partial_0^2 E^{(1)} \partial^k \partial^l E^{(1)} \\ &+ \left(3 + \left(c_s^{(2)}\right)^2\right) \partial_k \partial_l \partial_0 E^{(1)} \partial^k \partial^l \partial_0 E^{(1)} - \left(1 + \left(c_s^{(2)}\right)^2\right) \partial_m \partial_l \partial_k E^{(1)} \partial^m \partial^l \partial^k E^{(1)} - 2 \partial_0^2 \Psi^{(1)} \Delta E^{(1)}\\ &- \left(1 - 2 \left(c_s^{(2)}\right)^2\right) \partial_0 \Psi^{(1)} \Delta \partial_0 E^{(1)} + \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \partial_k \Delta E^{(1)} \partial^k \Psi^{(1)} + 2 \left(1 + \left(c_s^{(2)}\right)^2\right) \Delta E^{(1)} \Delta \Psi^{(1)}\\ &+ 2 \left(c_s^{(2)}\right)^2 \partial_k \partial_l \Psi^{(1)} \partial^k \partial^l E^{(1)} - 4{\cal{H}} \big( \partial_0 \Psi^{(1)} \left( 1 + \left(c_s^{(2)}\right)^2 \right) \Delta E^{(1)} \\ &+ \Psi^{(1)} \left(\left(c_s^{(2)}\right)^2 - c_s^2\right) \Delta \partial_0 E^{(1)} \big) \Biggr) \frac{1}{4} {\cal{T}}^{i j} \Biggr( -\Psi^{(1)} \partial_i \partial_j \Psi^{(1)} + \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_0 \Psi^{(1)} \partial_i \partial_j \partial_0 \Psi^{(1)}\\ &+\left(- 5 \partial^m \Psi^{(1)} \partial_m - 6 c_s^2 {\cal{H}} \partial_0 \Psi^{(1)} - 2 \left(1 - c_s^2\right) \Delta \Psi^{(1)} + \partial_0 \Psi^{(1)} \partial_0 - 2 \Psi^{(1)} \Delta \right) \partial_i \partial_j E^{(1)}+ 2 \partial_m \partial_0 E^{(1)} \partial^m \partial_i \partial_j \partial_0 E^{(1)}\\ & + \Delta \partial_0 E^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - \partial_m \Delta E^{(1)} \partial^m \partial_i \partial_j E^{(1)}- \partial_n \partial_m E^{(1)} \partial^n \partial^m \partial_i \partial_j E^{(1)} + 2 c_s^2 {\cal{H}} \Delta \partial_0 E^{(1)} \partial_i \partial_j E^{(1)} \Biggr) , \end{aligned} $

      (89)

      $ \begin{aligned}[b] \frac{1}{2} \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j} =& \frac{1}{2} \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) \Biggr( 3 \partial_i \Psi^{(1)} \partial_j \Psi^{(1)} - \frac{4}{3 (1 + w) {\cal{H}}^2} \partial_i \partial_0 \Psi^{(1)} \partial_j \partial_0 \Psi^{(1)} \\ &+ 2 \Psi^{(1)} \partial_i \partial_j \Psi^{(1)}- 2 \partial_k \partial_j \partial_0 E^{(1)} \partial^k \partial_i \partial_0 E^{(1)} - \partial_k \Delta E^{(1)} \partial^k \partial_i \partial_j E^{(1)} \\ &+ \partial_k \partial_l \partial_j E^{(1)} \partial^k \partial^l \partial_i E^{(1)} - 2 \Delta \partial_0^2 E^{(1)} \partial_i \partial_j E^{(1)} - 4{\cal{H}} \left(1 + c_s^2\right) \Delta \partial_0 E^{(1)} \partial_i \partial_j E^{(1)} \\ & + \Delta \partial_0 E^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - \partial_k \partial_i \partial_j E^{(1)} \partial^k \Psi^{(1)}+ 2 \partial_k \partial_j \Psi^{(1)} \partial^k \partial_i E^{(1)} + 2 \partial_k \partial_i \Psi^{(1)} \partial^k \partial_j E^{(1)} \\ &+ 12 \left(1 + c_s^2\right) {\cal{H}} \partial_0 \Psi^{(1)} \partial_i \partial_j E^{(1)} + 6 \partial_0^2 \Psi^{(1)} \partial_i \partial_j E^{(1)} \\ &- 4 (1 + c_s^2) \Delta \Psi^{(1)} \partial_i \partial_j E^{(1)} + \partial_0 \Psi^{(1)} \partial_i \partial_j \partial_0 E^{(1)} - 2 \Delta E^{(1)} \partial_i \partial_j \Psi^{(1)}\Biggr) . \end{aligned} $

      (90)

      It is not surprising that the left hand sides of the equations of motion are the same as those of Eqs. (68a)-(68d). As suggested in Ref. [47], we demonstrate that the source terms in Eq. (86b) include the first order perturbation $ E^{(1)} $ without a temporal derivative. This differs from the formulae in the other previous works [13, 44, 48].

    • 2.   Gauge invariant second order synchronous variables
    • For the gauge invariant synchronous metric perturbations in the second order, we also have

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} \mathrm{d} x^{\mu} \mathrm{d} x^{\nu} = & a^2 \left(- 2 \delta_{ij} \Psi^{(2)} + 2 \partial_i \partial_j E^{(2)} \right.\\ &\left.+ \partial_i C_j^{(2)} + \partial_j C_i^{(2)} + H_{ij}^{(2)}\right) \mathrm{d} x^i \mathrm{d} x^j. \end{aligned} $

      (91)

      The explicit expressions of the gauge invariant second order variables, $ E^{(2)} $, $ \Psi^{(2)} $, $ C_i^{(2)} $, and $ H_{ij}^{(2)} $ are presented in Eqs. (46a)-(46g), with the expression of $ {\cal{X}}_{\mu \nu}^S $ in Appendix E. Based on the temporal components of the second order Einstein field equations and Eqs. (59a), (59b), and (79a)-(80), the gauge invariant second order density perturbation is obtained as

      $ \begin{aligned}[b] \rho^{(\rm{GI}, 2)} = &- \frac{1}{3 (1 + w) \kappa a^2 {\cal{H}}^2} \Bigr(8 \partial_i \partial_0 \Psi^{(1)} \partial^i \partial_0 \Psi^{(1)} \\ &+ 6 (1 + w) {\cal{H}}^3 \Bigr(3 \partial_0 \Psi^{(2)} - \Delta \partial_0 E^{(2)} - 4 \partial_0 \Psi^{(1)} \Delta E^{(1)} \\ &+ 4 \Psi^{(1)} \Bigr(3 \partial_0 \Psi^{(1)} - \Delta \partial_0 E^{(1)}\Bigr) + 4 \partial_i \partial_j \partial_0 E^{(1)} \partial^i \partial^j E^{(1)}\Bigr) \\ &- 3 (1 + w) {\cal{H}}^2 \Bigr(6 \Bigr(\partial_0 \Psi^{(1)}\Bigr)^2 - 4 \partial_0 \Psi^{(1)} \Delta \partial_0 E^{(1)} \\ &+ 16 \Psi^{(1)} \Delta \Psi^{(1)} + 2 \Delta \Psi^{(2)} + 6 \partial_i \Psi^{(1)} \partial^i \Psi^{(1)} \\ &- 4 \partial_i \Delta E^{(1)} \partial^i \Psi^{(1)} + \Delta \partial_0 E^{(1)} \Delta \partial_0 E^{(1)} \\ &- 4 \Delta E^{(1)} \Delta \Psi^{(1)} - \partial_i \Delta E^{(1)} \partial^i \Delta E^{(1)} - 4 \partial_i \partial_j \Psi^{(1)} \partial^i \partial^j E^{(1)} \\ &- \partial_i \partial_j \partial_0 E^{(1)} \partial^i \partial^j \partial_0 E^{(1)} + \partial_i \partial_j \partial_k E^{(1)} \partial^i \partial^j \partial^k E^{(1)}\Bigr)\Bigr) . \end{aligned} $

      (92)

      Substituting Eqs. (59a), (59b), (79a)-(80), and (92) into the spatial parts of the gauge invariant second order Einstein field equations, we also obtain

      $ \begin{array}{l} {\cal{G}}_{i j} +{\cal{S}}_{i j} = 0, \end{array} $

      (93)

      where $ {\cal{G}}_{i j} $ is expressed in Eq. (66) and $ {\cal{S}}_{i j} $ in Eq. (85). Therefore, we derive the equations of motion of the gauge invariant second order cosmological perturbations as

      $ \partial_0^2 H_{i j}^{(2)} + 2{\cal{H}} \partial_0 H_{i j}^{(2)} - \Delta H_{i j}^{(2)} = - 4 \Lambda^{kl}_{ij} {\cal{S}}_{kl}, \tag{94a} $

      $ (\partial_0 + 2{\cal{H}}) \partial_0 C_i^{(2)} = - 4 \Delta^{- 1} {\cal{T}}^k_i \partial^l {\cal{S}}_{k l}, \tag{94b} $

      $ \left(\partial_0 + 2{\cal{H}}\right) \partial_0 E^{(2)} + \Psi^{(2)} = - 2 \Delta^{- 1} \left( \partial^i \Delta^{- 1} \partial^j - \frac{1}{2} {\cal{T}}^{i j} \right) {\cal{S}}_{i j}, \tag{94c} $

      $ \begin{aligned}[b] \partial_0^2 \Psi^{(2)} &+ \left(2 + 3 \left(c_s^{(2)}\right)^2\right) {\cal{H}} \partial_0 \Psi^{(2)} - \frac{1}{2} \left(1 + 2 \left(c_s^{(2)}\right)^2\right) \Delta \Psi^{(2)} \\ &- \Delta \left( \frac{1}{2} \left(\partial_0 + 2{\cal{H}}\right) \partial_0 E^{(2)} + \left(c_s^{(2)}\right)^2 {\cal{H}} \partial_0 E^{(2)} \right) \\ &= - \frac{1}{2} {\cal{T}}^{i j} {\cal{S}}_{i j} . \end{aligned} \tag{94d}$

      The right hand sides of the above equations are the same as that in Eqs. (87)-(90), because the source term $ {\cal{S}}_{i j} $ is uniquely determined by the gauge invariant first order synchronous variables. In particular, the equation of motion of the tensor perturbation $ H_{\mu \nu}^{(\rm{GI},2)} $ in Eq. (94a) is the same as that in Eq. (86a). For the vector perturbation, the evolution of $ \partial_0 C_i^{(2)} $ in Eq. (94b) is demonstrated to be as same as that of $ V_i^{(2)} $ in Eq. (86b).

    VI.   CONCLUSIONS AND DISCUSSIONS
    • In this study, we investigated the gauge invariance of the cosmological perturbations up to the second order via the Lie derivative method. We demonstrated that there are infinite families of gauge invariant variables for the cosmological perturbations up to the second order. For different families, we determined their conversion formulae, which belong to a linear space spanned by a finite number of bases that were also shown to be gauge invariant. In particular, we focused on the Newtonian, synchronous, and hybrid variables. We presented the explicit conversions between these different families of gauge-invariant variables. In contrast to the first order, the second order gravitational waves were demonstrated to be mixed with the first order cosmological perturbations. Therefore, the gauge invariance is important in the studies on second order gravitational waves.

      We derived the equations of motion of the gauge invariant second order cosmological perturbations, which were obtained from the gauge invariant first order scalar perturbations. It was determined that the choices of gauge-invariant variables at different orders are independent,

      e.g., Newtonian at first order while synchronous at second order, and vice versa. In this research, we studied both of the above four typical cases. To obtain the gauge invariant equations of motion, we decomposed the gauge invariant perturbed Einstein field equations into scalar, vector, and tensor components.

      In principle, the above method can be generalized to explore the higher order cosmological perturbations. A formal derivation of the gauge invariant higher order cosmological perturbations has been demonstrated in a previous study [18]. Based on the same approach in this work, it is straightforward to obtain the explicit expressions for the gauge invariant higher order scalar, vector, and tensor perturbations, as well as their equations of motion.

      This study was based on the method developed by [5, 7]. We further focused on the conversion between different gauge-invariant variables. Because there is no uniqueness in the framework, it might suggest that one can not determine the gauge-invariant tensor perturbations $ H^{(2)}_{i j} $ that correspond to the energy density spectrum of gravitational waves. However, physical insights or arguments should be explored in future studies.

    ACKNOWLEDGMENTS
    • We acknowledge Prof. Rong-Gen Cai, Prof. Qing-Guo Huang, Prof. Xin Li, Prof. Tao Liu, Prof. Shi Pi, Prof. Jiang-Hao Yu, Prof. Xin Zhang, Mr. Zu-Cheng Chen, Mr. Chen Yuan, Mr. Xukun Zhang, and Mr. Jing-Zhi Zhou for their helpful discussions. We acknowledge the xPand package [59].

    APPENDIX A: EXPANSION OF A GENERIC TENSOR WITH LIE DERIVATIVE
    • For a generic tensor $ {\cal{Q}} $, its Lie derivative along a vector $ \zeta^{\mu} $ is defined as

      $ {\cal{L}}_{\zeta} {\cal{Q}} \equiv \lim_{\epsilon \rightarrow 0} \frac{\varphi^{\ast}_{\epsilon} {\cal{Q}} (x) -{\cal{Q}} (x)}{\epsilon}, \tag{A1}$

      where $ \varphi^{\ast}_{\epsilon} {\cal{Q}} $ is a one-parameter coordinate transformation that acts on the tensor $ {\cal{Q}} $. In the $ \epsilon \rightarrow 0 $ regime, it becomes the infinitesimal transformation in the form of Eq. (1), where $ \xi^{(1) \mu} \propto \epsilon \zeta^{\mu} $. For a scalar function, contravariant vector, and covariant vector, the expressions of $ \varphi_{\epsilon}^{\ast} $ are given by

      $ \varphi^{\ast}_{\epsilon} f (x) = f (\tilde{x}), \tag{A2a} $

      $ \varphi^{\ast}_{\epsilon} w_{\mu} (x) = \frac{\partial \tilde{x}^{\nu}}{\partial x^{\mu}} w_{\nu} (\tilde{x}), \tag{A2b} $

      $ \varphi_{\epsilon, \ast} A^{\mu} (x) = \frac{\partial x^{\mu}}{\partial \tilde{x}^{\nu}} A^{\nu} (\tilde{x}) .\tag{A2c} $

      Based on the first order expansions of Eqs. (A2a)-(A2c), we obtain

      $ {\cal{L}}_{\zeta} f = \zeta^{\nu} \partial_{\nu} f, \tag{A3a} $

      $ {\cal{L}}_{\zeta} w_{\mu} = \zeta^{\nu} \partial_{\nu} w_{\mu} + w_{\nu} \partial_{\mu} \zeta^{\mu},\tag{A3b} $

      $ {\cal{L}}_{\zeta} A^{\mu} = \zeta^{\nu} \partial_{\nu} A^{\mu} - A^{\nu} \partial_{\nu} \zeta^{\mu} . \tag{A3c} $

      The Lie derivative in Eq. (A3c) is also symbolized as $ {\cal{L}}_X A^{\mu} \equiv [X, A]^{\mu} $, where $ [,] $ is the Lie bracket. For a tensor $ S_{\mu \nu} $, its Lie derivative is given by

      $ {\cal{L}}_{\zeta} S_{\mu \nu} = \zeta^{\lambda} \partial_{\lambda} S_{\mu \nu} + S_{\lambda \nu} \partial_{\mu} \zeta^{\lambda} + S_{\mu \lambda} \partial_{\nu} \zeta^{\lambda} . \tag{A4} $

      Higher order expansions of any tensor upon the infinitesimal transformation have been constructed in terms of Lie derivatives [5-7, 10, 11]. In the following, we briefly review the formulae up to the second order.

      For the scalar function $ f (x) $ upon the infinitesimal transformation, up to the second order, we expand it in terms of $ \xi^{(1)} $ and $ \xi^{(2)} $, namely,

      $ \begin{aligned}[b] f (x) \rightarrow &f (\tilde{x}) \approx f \left( x + \xi^{(1)} + \frac{1}{2} \xi^{(2)} \right) \\ = &f + \xi^{(1) \nu} \partial_{\nu} f + \frac{1}{2} \xi^{(2) \nu} \partial_{\nu} f + \frac{1}{2} \xi^{(1) \sigma} \xi^{(1) \rho} \partial_{\sigma} \partial_{\rho} f \\ = &f + \xi^{(1) \nu} \partial_{\nu} f + \frac{1}{2} \left.(\xi^{(2) \nu} \partial_{\nu} f + \xi^{(1) \sigma} \xi^{(1) \rho} \partial_{\sigma} \partial_{\rho} f \right.\\ &\left.+ \xi^{(1) \sigma} \partial_{\sigma} \xi^{(1) \rho} \partial_{\rho} f - \xi^{(1) \sigma} \partial_{\sigma} \xi^{(1) \rho} \partial_{\rho} f\right) \\ = & f +{\cal{L}}_{\xi^{(1)}} f + \frac{1}{2} \left({\cal{L}}_{\xi^{(2)} - \xi^{(1) \nu} \partial_{\nu} \xi^{(1)}} +{\cal{L}}^2_{\xi^{(1)}}\right) f. \end{aligned} \tag{A5} $

      For simplicity, we denote

      $ \xi_1^{\mu} \equiv \xi^{(1) \mu}, \tag{A6a}$

      $ \xi_2^{\mu} \equiv \xi^{(2) \mu} - \xi^{(1) \nu} \partial_{\nu} \xi^{(1) \mu} . \tag{A6b} $

      It leads to

      $ f (\tilde{x}) = f +{\cal{L}}_{\xi_1} f + \frac{1}{2} \left({\cal{L}}_{\xi_2} +{\cal{L}}^2_{\xi_1}\right) f +{\cal{O}} \left(\xi^{(3)}\right) . \tag{A7}$

      For the contravariant vector $ w_{\mu} $, up to the second order, we expand it in terms of $ \xi^{(1)} $ and $ \xi^{(2)} $, namely,

      $ \begin{aligned}[b] w_{\mu} (x) \rightarrow & \frac{\partial \tilde{x}^{\nu}}{\partial x^{\mu}} w_{\nu} (\tilde{x}) \approx \left( \delta^{\lambda}_{\mu} + \partial_{\mu} \xi^{(1) \lambda} + \frac{1}{2} \partial_{\mu} \xi^{(2) \lambda} \right) \Biggr( w_{\lambda} (x) \\ & + \left( \xi^{(1) \nu} + \frac{1}{2} \xi^{(2) \nu} \right) \partial_{\nu} w_{\lambda} \\ &+ \left( \xi^{(1) \sigma} + \frac{1}{2} \xi^{(2) \sigma} \right) \left( \xi^{(1) \rho} + \frac{1}{2} \xi^{(1) \rho} \right) \partial_{\sigma} \partial_{\rho} w_{\lambda} \Biggr) \\ =& w_{\mu} (x) + \xi^{(1) \nu} \partial_{\nu} w_{\mu} + w_{\lambda} \partial_{\mu} \xi^{(1) \lambda} + \frac{1}{2} w_{\lambda} \partial_{\mu} \xi^{(2) \lambda} \\ &+ \partial_{\mu} \xi^{(1) \lambda} \xi^{(1) \nu} \partial_{\nu} w_{\lambda} \\ &+ \frac{1}{2} \left(\left(\xi^{(2) \nu} \partial_{\nu} w_{\mu} + \xi^{(1) \sigma} \xi^{(1) \rho} \partial_{\sigma} \partial_{\rho} w_{\mu}\right)\right) \\ = & w_{\mu} +{\cal{L}}_{\xi^{(1)}} w_{\mu} + \frac{1}{2} {\cal{L}}_{\xi^{(2)}} w_{\mu} \\ &+ \frac{1}{2} \left({\cal{L}}_{\xi^{(1)}}^2 w_{\mu} - \xi^{(1) \sigma} \partial_{\sigma} \xi^{(1) \rho} \partial_{\rho} w_{\mu} \right)\\ &- w_{\rho} \partial_{\mu} \left(\left(\xi^{(1) \sigma} \partial_{\sigma} \xi^{(1) \rho}\right)\right) \\ =& w_{\mu} +{\cal{L}}_{\xi^{(1)}} w_{\mu} + \frac{1}{2} \left({\cal{L}}_{\xi^{(2)} - \xi^{(1) \sigma} \partial_{\sigma} \xi^{(1)}} +{\cal{L}}_{\xi^{(1)}}^2\right) w_{\mu} \\ =& w_{\mu} +{\cal{L}}_{\xi_1} w_{\mu} + \frac{1}{2} \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right) w_{\mu} . \end{aligned} \tag{A8}$

      Similarly, for the covariant vector $ A^{\mu} $, up to the second order, we expand it in terms of $ \xi^{(1)} $ and $ \xi^{(2)} $, namely,

      $ \begin{aligned}[b] A^{\mu} (x) \rightarrow & \frac{\partial x^{\mu}}{\partial \tilde{x}^{\nu}} A^{\nu} (\tilde{x}) \approx \left( \frac{\partial \tilde{x}^{\mu}}{\partial x^{\nu}} \right)^{- 1} A^{\nu} \left( x + \xi^{(1)} + \frac{1}{2} \xi^{(2)} \right) \\ \approx & A^{\mu} + \xi^{(1) \sigma} \partial_{\sigma} A^{\mu} - A^{\nu} \partial_{\nu} \xi^{(1) \mu} - \partial_{\nu} \xi^{(1) \mu} \xi^{(1) \sigma} \partial_{\sigma} A^{\nu} \\ &- \frac{1}{2} \left(\partial_{\nu} \xi^{(2) \mu} - 2 \partial_{\lambda} \xi^{(1) \mu} \partial_{\nu} \xi^{(1) \lambda}\right) A^{\nu} \\ &+ \frac{1}{2} \Biggr(\left(\xi^{(2) \sigma} \partial_{\sigma} A^{\mu} + \xi^{(1) \sigma} \xi^{(1) \rho} \partial_{\sigma} \partial_{\rho} A^{\mu}\right) \\ =& A^{\mu} +{\cal{L}}_{\xi^{(1)}} A^{\mu} + \frac{1}{2} \left({\cal{L}}_{\xi^{(2)} - \xi^{(1) \nu} \partial_{\nu} \xi^{(1)})} A^{\mu} +{\cal{L}}^2_{\xi^{(1)}} A^{\mu}\right)\Biggr) \\ =& A^{\mu} +{\cal{L}}_{\xi_1} A^{\mu} + \frac{1}{2} \left({\cal{L}}_{\xi_2} A^{\mu} +{\cal{L}}^2_{\xi_1} A^{\mu}\right), \end{aligned} \tag{A9} $

      where the inverse Jacobi matrix $ {\partial x^{\mu}}/{\partial \tilde{x}^{\nu}} $ is given by

      $ \begin{aligned}[b] \left( \frac{\partial \tilde{x}^{\mu}}{\partial x^{\nu}} \right)^{- 1} =& \delta^{\mu}_{\nu} - \partial_{\nu} \xi^{(1) \mu} - \frac{1}{2} \left(\partial_{\nu} \xi^{(2) \mu}\right.\\ &\left.- 2 \partial_{\lambda} \xi^{(1) \mu} \partial_{\nu} \xi^{(1) \lambda}\right) +{\cal{O}} \left(\xi^{(3)}\right) . \end{aligned}\tag{A10} $

      It is derived from $ \left( \dfrac{\partial \tilde{x}^{\lambda}}{\partial x^{\nu}} \right)^{- 1} \dfrac{\partial \tilde{x}^{\mu}}{\partial x^{\lambda}} = \dfrac{\partial x^{\lambda}}{\partial \tilde{x}^{\nu}} \dfrac{\partial \tilde{x}^{\mu}}{\partial x^{\lambda}} = \delta^{\mu}_{\nu} $. Finally, we consider the tensor $ S_{\mu \nu} $, which is expanded as

      $ \begin{aligned}[b] S_{\sigma \rho} (x) \rightarrow & \frac{\partial \tilde{x}^{\mu}}{\partial x^{\sigma}} \frac{\partial \tilde{x}^{\nu}}{\partial x^{\rho}} S_{\mu \nu} (\tilde{x}) \approx S_{\mu \nu} \left( x + \xi^{(1)} + \frac{1}{2} \xi^{(2)} \right) \partial_{\sigma} \Biggr( x^{\mu} \\ &+ \xi^{(1) \mu} + \frac{1}{2} \xi^{(2) \mu} \Biggr) \partial_{\rho} \left( x^{\nu} + \xi^{(1) \nu} + \frac{1}{2} \xi^{(2) \nu} \right) \\ \approx & S_{\sigma \rho} + \xi^{(1) \lambda} \partial_{\lambda} S_{\sigma \rho} + S_{\sigma \nu} \partial_{\rho} \xi^{(1) \nu} + S_{\mu \rho} \partial_{\sigma} \xi^{(1) \mu}\\ &+ \frac{1}{2} \left(\left(\xi^{(2) \lambda} \partial_{\lambda} S_{\sigma \rho} + S_{\mu \rho} \partial_{\sigma} \xi^{(2) \mu} + S_{\sigma \nu} \partial_{\rho} \xi^{(2) \nu} \right.\right.\\ &+ \xi^{(1) \lambda} \xi^{(1) \kappa} \partial_{\lambda} \partial_{\kappa} S_{\sigma \rho} + 2 S_{\mu \nu} \partial_{\sigma} \xi^{(1) \mu} \partial_{\rho} \xi^{(1) \nu} \\ &\left.+ 2 \xi^{(1) \lambda} \partial_{\lambda} S_{\sigma \nu} \partial_{\rho} \xi^{(1) \nu} + 2 \xi^{(1) \lambda} \partial_{\lambda} S_{\mu \rho} \partial_{\sigma} \xi^{(1) \mu}\right) \\ = &S_{\sigma \rho} +{\cal{L}}_{\xi^{(1)}} S_{\sigma \rho} + \frac{1}{2} \left({\cal{L}}_{\xi^{(2)} - \xi^{(1) \nu} \partial_{\nu} \xi^{(1)}} +{\cal{L}}_{\xi^{(1)}}^2\right) S_{\sigma \rho} \\ =& S_{\sigma \rho} +{\cal{L}}_{\xi_1} S_{\sigma \rho} + \frac{1}{2} \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right) S_{\sigma \rho} . \end{aligned} \tag{A11}$

      The above formulae can be straightforwardly generalized to higher order expansions.

      In these derivations, we adopt Eq. (1) to expand $ \tilde{x} $. It can be obtained via

      $ \tilde{x}^\mu \equiv \varphi_\epsilon x^\mu = x^\mu + \frac{{\rm d}\varphi_\epsilon x^\mu}{{\rm d} \epsilon}\left|_{\epsilon = 0}\epsilon + \frac{1}{2}\frac{{\rm d}^2\varphi_\epsilon x^\mu}{{\rm d} \epsilon^2}\right|_{\epsilon = 0}\epsilon^2 + {\cal{O}}\left(\epsilon^3\right)\; , \tag{A12}$

      where the explicit expressions of $ \xi^{(1)} $ and $ \xi^{(2)} $ are defined as

      $ \xi^{(1)}\equiv\frac{{\rm d}\varphi_\epsilon x^\mu}{{\rm d} \epsilon}\Biggr|_{\epsilon = 0}\epsilon\; , \tag{A13a}$

      $\xi^{(2)}\equiv \frac{{\rm d}^2\varphi_\epsilon x^\mu}{{\rm d} \epsilon^2}\Biggr|_{\epsilon = 0}\epsilon^2\; .\tag{A13b} $

      Therefore, the order of the $ \xi^{(n)} $ is determined by the power of $ \epsilon $ in the transformation (Eq. (A12)), and has no relevance with metric perturbations or matter perturbations in principle.

    APPENDIX B: GAUGE TRANSFORMATIONS IN THE LANGUAGE OF LIE DERIVATIVE
    • For an arbitrary tensor ${\cal{Q}} (x^{\mu})^{i_1 i_2 \cdots}_{j_1 j_2 \cdots}$ upon the infinitesimal transformation, it can be expanded as [5]

      $ \begin{aligned}[b] {\cal{Q}} (x)^{i_1 i_2 \cdots}_{j_1 j_2 \cdots} \rightarrow & \left( \frac{\partial \tilde{x}^{l_1}}{\partial {x}^{j_1}} \frac{\partial \tilde{x}^{l_2}}{\partial {x}^{j_2}} \frac{\partial {x}^{i_1}}{\partial \tilde{x}^{k_1}} \frac{\partial {x}^{i_2}}{\partial \tilde{x}^{k_2}} \cdots \right) {\cal{Q}} (\tilde{x})^{k_1 k_2 \cdots}_{l_1 l_2 \cdots} \approx \left( \delta^{l_1}_{j_1} + \partial_{j_1} \xi^{(1) l_1} + \frac{1}{2} \partial_{j_1} \xi^{(2) l_1} \right) \left(\cdots\right) \Biggr({\cal{Q}} (x)^{k_1 k_2 \cdots}_{l_1 l_2 \cdots} + \xi^{(1) \nu} \partial_{\nu} {\cal{Q}} (x)^{k_1 k_2 \cdots}_{l_1 l_2 \cdots} \\ & + \frac{1}{2} \left((\xi^{(1) \sigma} \xi^{(2) \rho} \partial_{\sigma} \partial_{\rho} {\cal{Q}}(x)^{k_1 k_2 \cdots}_{l_1 l_2 \cdots} + \xi^{(2) \nu} \partial_{\nu} {\cal{Q}}(x)^{k_1 k_2 \cdots}_{l_1 l_2 \cdots}) \right) \Biggr) \approx {\cal{Q}} (x^{\mu})^{i_1 i_2 \cdots}_{j_1 j_2 \cdots} +{\cal{L}}_{\xi^{(1)}} {\cal{Q}} (x^{\mu})^{i_1 i_2 \cdots}_{j_1 j_2 \cdots} + \frac{1}{2} ({\cal{L}}_{\xi^{(2)} - \xi^{(1) \nu} \partial_{\nu} \xi^{(1)}} +{\cal{L}}_{\xi^{(1)}}^2) {\cal{Q}} (x^{\mu})^{i_1 i_2 \cdots}_{j_1 j_2 \cdots} \end{aligned} \tag{B1}$

      For simplicity, the above equation can be rewritten without indices as follows

      $ {\cal{Q}} \rightarrow \tilde{{\cal{Q}}} = {\cal{Q}}+{\cal{L}}_{\xi_1} {\cal{Q}}+ \frac{1}{2} ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) {\cal{Q}}+{\cal{O}} (\xi^{(3)}) \tag{B2} $

      In Appendix A, we presented four examples of Eq. (B2): Eqs. (A5), (A8), (A9), and (A11).

      In the following, we consider the infinitesimal transformation of perturbations up to the second order. The n-th order perturbations $ {\cal{Q}}^{(n)} $ of the tensor $ {\cal{Q}} $ can be expanded as follows

      $ {\cal{Q}} = {\cal{Q}}^{(0)} +{\cal{Q}}^{(1)} + \frac{1}{2} {\cal{Q}}^{(2)} +{\cal{O}} ({\cal{Q}}^{(3)}) . \tag{B3} $

      In order to study the gauge transformation of $ {\cal{Q}}^{(n)} $, one can choose a type of infinitesimal transformation, where the order of $ \xi^{(1)} $ is in the same order as $ {\cal{Q}}^{(1)} $. Combining Eqs. (B2) with (B3), we obtain

      $ \begin{aligned}[b] \tilde{{\cal{Q}}} =& {\cal{Q}}^{(0)} + \left({\cal{Q}}^{(1)} +{\cal{L}}_{\xi_1} {\cal{Q}}^{(0)}\right) \\ &+ \frac{1}{2} \left({\cal{Q}}^{(2)} + 2{\cal{L}}_{\xi_1} {\cal{Q}}^{(1)} + \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right){\cal{Q}}^{(0)}\right) +{\cal{O}} \left(\xi^{(3)}\right) . \end{aligned} \tag{B4}$

      On the other side, the tensor $ \tilde{{\cal{Q}}} $ can also be expanded in the form of Eq. (B3), namely,

      $ \tilde{{\cal{Q}}} = \tilde{{\cal{Q}}}^{(0)} + \tilde{{\cal{Q}}}^{(1)} + \frac{1}{2} \tilde{{\cal{Q}}}^{(2)} +{\cal{O}} (\tilde{{\cal{Q}}}^{(3)} ) . \tag{B5} $

      Therefore, based on Eqs. (B4) and (B5), we deduce the infinitesimal transformations of the zeroth, first, and second order perturbations in Eqs. (2).

    APPENDIX C: BRIEF DERIVATION OF EQS. (5A) AND (5B)
    • For the first order metric perturbations $ g^{(1)}_{\mu \nu} $, it could be split into a gauge invariant part $ g^{(\mathrm{\rm{GI}, 1)}}_{\mu \nu} $ and a gauge variant counter term $ \mathbb{C}_{\mu \nu}^{(1)} $, i.e.,

      $ g_{\mu \nu}^{(1)} = : g_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} + \mathbb{C}_{\mu \nu}^{(1)} . \tag{C1} $

      Based on the gauge transformation of $ g_{\mu \nu}^{(1)} $ in Eq. (3a), we could rewrite $ g^{(\mathrm{\rm{GI}, 1)}}_{\mu \nu} $ to be

      $ \begin{aligned}[b] g^{(\mathrm{\rm{GI}, 1)}}_{\mu \nu} = \tilde{g}_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}}= \tilde{g}_{\mu \nu}^{(1)} - \tilde{\mathbb{C}}_{\mu \nu}^{(1)} = g_{\mu \nu}^{(1)} +{\cal{L}}_{\xi_1} g_{\mu \nu}^{(0)} - \tilde{\mathbb{C}}_{\mu \nu}^{(1)} . \end{aligned}\tag{C2} $

      Based on Eqs. (C1) and (C2), we obtain

      $ \tilde{\mathbb{C}}^{(1)}_{\mu \nu} = \mathbb{C}^{(1)}_{\mu \nu} +{\cal{L}}_{\xi_1} g_{\mu \nu}^{(0)} . \tag{C3} $

      If we let $ \mathbb{C}^{(1)}_{\mu \nu} \equiv {\cal{L}}_X g_{\mu \nu}^{(0)} $, Eq. (C3) is reduced to ${\cal{L}}_{\tilde{X}} g_{\mu \nu}^{(0)} = $$ {\cal{L}}_{X + \xi_1} g_{\mu \nu}^{(0)}$. We can obtain

      $ \tilde{X}^{\mu} - X^{\mu} = \xi_1^{\mu} . \tag{C4} $

      Here, $ X^\mu $ is rewritten in a form independent of the background metric. Therefore, we can rewrite Eq. (C1) as

      $ g_{\mu \nu}^{(1)} = g_{\mu \nu}^{(\mathrm{\rm{GI}, 1)}} +{\cal{L}}_X g_{\mu \nu}^{(0)} . \tag{C5} $

      For the second order metric perturbations $ g_{\mu \nu}^{(2)} $, we could also split it as

      $ g_{\mu \nu}^{(2)} = : g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} + \mathbb{C}_{\mu \nu}^{(2)}, \tag{C6} $

      where $ g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} $ are the gauge invariant second order metric perturbations. Based on Eqs. (3b) and (C6), we can also express $ g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}} $ as

      $ \begin{aligned}[b] g_{\mu \nu}^{(\mathrm{\rm{GI}, 2)}}& = \tilde{g}_{\mu \nu}^{(\mathrm{\rm{GI},2)}} = \tilde{g}_{\mu \nu}^{(2)} - \tilde{\mathbb{C}}_{\mu \nu}^{(2)} \\ &= 2{\cal{L}}_{\xi_1} g_{\mu \nu}^{(1)} + ({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2) g_{\mu \nu}^{(0)} - \tilde{\mathbb{C}}_{\mu \nu}^{(2)} . \end{aligned} \tag{C7}$

      Based on Eqs. (C6) and (C7), we have

      $ \begin{aligned}[b] \tilde{\mathbb{C}}_{\mu \nu}^{(2)} - \mathbb{C}_{\mu \nu}^{(2)} =& 2{\cal{L}}_{\xi_1} g_{\mu \nu}^{(1)} + \left({\cal{L}}_{\xi_2} +{\cal{L}}_{\xi_1}^2\right) g_{\mu \nu}^{(0)} = 2{\cal{L}}_{\tilde{X}} \left(\tilde{g}_{\mu \nu}^{(1)} -{\cal{L}}_{\xi_1} g_{\mu \nu}^{(0)}\right) - 2{\cal{L}}_X g_{\mu \nu}^{(1)} +{\cal{L}}_{\xi_2} g_{\mu \nu}^{(0)} + \left({\cal{L}}_{\tilde{X}}^2 +{\cal{L}}_X^2 -{\cal{L}}_{\tilde{X}} {\cal{L}}_X -{\cal{L}}_X {\cal{L}}_{\tilde{X}}\right) g_{\mu \nu}^{(0)} \\ = &\left(2{\cal{L}}_{\tilde{X}} \tilde{g}_{\mu \nu}^{(1)} - {\cal{L}}^2_{\tilde{X}} g_{\mu \nu}^{(0)}\right) - \left(2{\cal{L}}_X g_{\mu \nu}^{(1)} -{\cal{L}}_X^2 g_{\mu \nu}^{(0)}\right) +{\cal{L}}_{\left(\xi_2 + [\tilde{X}, X]\right)} g_{\mu \nu}^{(0)} \equiv \left(2{\cal{L}}_{\tilde{X}} \tilde{g}_{\mu \nu}^{(1)} + \left({\cal{L}}_{\tilde{Y}} - {\cal{L}}^2_{\tilde{X}}\right) g_{\mu \nu}^{(0)}\right) \\&- \left(2{\cal{L}}_X g_{\mu \nu}^{(1)} + \left({\cal{L}}_Y -{\cal{L}}_X^2\right) g_{\mu \nu}^{(0)}\right), \end{aligned} \tag{C8}$

      where we have introduced the infinitesimal vector $ Y^{\mu} $ satisfying

      $ \tilde{Y}^{\mu} - Y^{\mu} = \xi_2^{\mu} + [\xi_1, X]^{\mu} . \tag{C9} $

      Therefore, we can rewrite Eq. (C6) in the form of

      $ g_{\mu \nu}^{(2)} = g^{(\mathrm{\rm{GI}, 2)}}_{\mu \nu} + 2{\cal{L}}_X g^{(1)}_{\mu \nu} + ({\cal{L}}_Y -{\cal{L}}_X^2) g_{\mu \nu}^{(0)} .\tag{C10} $

    APPENDIX D: SCALAR-VECTOR-TENSOR DECOMPOSITION
    • By adopting the transverse operators, namely, $ {\cal{T}}^i_j \equiv \delta^i_j - \partial^i \Delta^{- 1} \partial_j $, we could decompose an arbitrary spatial vector $ U_i $ to be

      $ U_i = U^{(T)}_i + \partial_i U^{(S)}, \tag{D1} $

      where the transverse $ U_i^{(T)} $ and longitudinal $ U^{(S)} $ parts are defined as

      $ U_i^{(T)} \equiv {\cal{T}}^k_i U_k, \tag{D2a} $

      $ U^{(S)} \equiv \Delta^{- 1} \partial^i U_i , \tag{D2b}$

      and $ \Delta^{- 1} $ is the inverse Laplacian operator.

      For an arbitrary spatial tensor $ S_{ij} $, we could decompose it as [60, 61]

      $ \begin{aligned}[b] S_{ij} =& S^{(H)}_{ij} + 2 \delta_{ij} S^{(\Psi)} + 2 \partial_i \partial_j S^{(E)} \\&+ \partial_j S_i^{(C)} + \partial_i S_j^{(C)},\end{aligned} \tag{D3} $

      where we can define the tensor $ S^{(H)}_{ij} $, vector $ S_i^{(C)} $, and scalar $ S^{(\Psi)} $, $ S^{(E)} $ parts as

      $ S^{(H)}_{ij} \equiv \left( {\cal{T}}^k_i {\cal{T}}^l_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{kl} \right) S_{kl}, \tag{D4a} $

      $ S_i^{(C)} \equiv \Delta^{- 1} \partial^l {\cal{T}}^k_i S_{lk} , \tag{D4b}$

      $ S^{(\Psi)} \equiv \frac{1}{4} {\cal{T}}^{kl} S_{kl}, \tag{D4c}$

      $ S^{(E)} \equiv \frac{1}{2} \Delta^{- 1} \left( \partial^k \Delta^{- 1} \partial^l - \frac{1}{2} {\cal{T}}^{kl} \right) S_{lk}. \tag{D4d}$

      (127d)

      Here, ($ {\cal{T}}^k_i {\cal{T}}^l_j - {\cal{T}}_{ij} {\cal{T}}^{kl}/2 $) denotes the transverse and traceless operator. Other kinds of decomposition can be found in Refs. [52, 62].

    APPENDIX E: EXPLICIT EXPRESSIONS OF $ {\cal{X}}_{\mu \nu} $
    • For a given $ X^{\mu} $, we express $ {\cal{X}}_{\mu \nu} $ in terms of the first order metric perturbations, based on Eq. (40). The obtained results are expressed as

      $ \begin{aligned}[b] {\cal{X}}_{\mu \nu} =& \frac{4}{a^2}{\cal{H}}X^0 {g}^{(1)}_{\mu \nu} + 2 X^{\sigma} \partial_{\sigma} \left(\frac{1}{a^2}{g}_{\mu \nu}^{(1)}\right) + \frac{2}{a^2} {g}_{\sigma \nu}^{(1)} \partial_{\mu} X^{\sigma} \\ &+ \frac{2}{a^2} {g}_{\mu \sigma}^{(1)} \partial_{\nu} X^{\sigma}- \eta_{\mu \nu} \Big(\Big(4{\cal{H}}^2 + 2 \dot{{\cal{H}}}\Big) \Big(X^0\Big)^2 \\ &+ 2{\cal{H}}X^{\mu} \partial_{\mu} X^0\Big) - 2 \eta_{\sigma \rho} \partial_{\mu} X^{\sigma} \partial_{\nu} X^{\rho} - \eta_{\mu \sigma} \Big(X^{\rho} \partial_{\rho} \partial_{\nu} X^{\sigma} \\ &+ \partial_{\rho} X^{\sigma} \partial_{\nu} X^{\rho} + 4{\cal{H}}X^0 \partial_{\nu} X^{\sigma}\Big) - \eta_{\sigma \nu} \Big(X^{\rho} \partial_{\rho} \partial_{\mu} X^{\sigma} \\ &+ \partial_{\rho} X^{\sigma} \partial_{\mu} X^{\rho} + 4{\cal{H}}X^0 \partial_{\mu} X^{\sigma}\Big) \end{aligned}\tag{E1} $

        For $ X^{\mu} $ in Eq. (22), we obtain $ {\cal{X}}_{\mu \nu} $ to be

      $ \begin{aligned}[b] \mathcal{X}_{00}^N = & - 2 \left( \frac{2}{a} (\partial_0 (a (\partial_0 {e}^{(1)} - b^{(1)}))) + (\partial_0 {e}^{(1)} - b^{(1)}) \partial_0 + \delta^{ik} (c_k^{(1)} + \partial_k {e}^{(1)}) \partial_i \right) \Biggr( 2 \phi^{(1)} - \frac{1}{a} \partial_0 \left(a (\partial_0 {e}^{(1)} - b^{(1)})\right) \Biggr) \\ & + 2 \delta^{ik} (\partial_0 \left(c_k^{(1)} + \partial_k {e}^{(1)}) \right) \left(\partial_i b^{(1)} + 2 \nu_i^{(1)} - \partial_0 c_i^{(1)}\right), \end{aligned} \tag{E2a}$

      $ \begin{aligned}[b] \mathcal{X}_{0 j}^N = & \Biggr( \delta^k_j \Biggr( \frac{2 \dot{a}}{a} \left(\partial_0 {e}^{(1)} - b^{(1)}\right) + \left(\partial_0 {e}^{(1)} - b^{(1)}\right) \partial_0 + \delta^{ik} \left(c_k^{(1)} + \partial_k {e}^{(1)}\right) \partial_i + \left(\partial_0 \left(\partial_0 {e}^{(1)} - b^{(1)}\right)\right) \Biggr) \\ &+ \delta^{ks} \left(\partial_j \left(c_s^{(1)} + \partial_s {e}^{(1)}\right)\right) \Biggr) \left(\partial_k b^{(1)} + 2 \nu_k^{(1)} - \partial_0 c_k^{(1)}\right) - 2 \left( \partial_j \left(\partial_0 {e}^{(1)} - b^{(1)}\right) \right) \Biggr( 2 \phi^{(1)} - \frac{1}{a} \partial_0 \left(a \left(\partial_0 {e}^{(1)} - b^{(1)}\right)\right) \Biggr) \\ & + \delta^{kl} \left(\partial_0 \left(c_l^{(1)} + \partial_l {e}^{(1)}\right)\right) \Biggr( - 2 \delta_{jk} \Biggr( 2 \psi^{(1)} + \frac{\dot{a}}{a} (\partial_0 {e}^{(1)} - b^{(1)}) \Biggr) + 2 \partial_j \partial_k {e}^{(1)} + \partial_j c_k^{(1)} + \partial_k c_j^{(1)} + 2h_{j k}^{(1)} \Biggr), \end{aligned} \tag{E2b}$

      (129b)

      $ \begin{aligned}[b] \mathcal{X}_{kl}^N = & \Biggr(\delta^s_k \delta^t_l \Biggr( \frac{2 \dot{a}}{a} \left(\partial_0 {e}^{(1)} - b^{(1)}\right) + \left(\partial_0 {e}^{(1)} - b^{(1)}\right) \partial_0 + \delta^{ik} \left(c_k^{(1)} + \partial_k {e}^{(1)}\right) \partial_i \Biggr) \\ &+ \left(\left(\delta^s_k \delta^{tw} \partial_l + \delta^t_l \delta^{sw} \partial_k\right) \left(c_w^{(1)} + \partial_w {e}^{(1)}\right)\right)\Biggr) \Biggr( - 2 \delta_{st} \Biggr( 2 \psi^{(1)} + \frac{\dot{a}}{a} \left(\partial_0 {e}^{(1)} - b^{(1)}\right) \Biggr) + 2 \partial_s \partial_t {e}^{(1)} \\ &+ \partial_s c_t^{(1)} + \partial_t c_s^{(1)} + 2h_{s t}^{(1)} \Biggr) + \left(\left(\delta^s_k \partial_l + \delta^s_l \partial_k) (\partial_0 {e}^{(1)} - b^{(1)}\right)\right) \left(\partial_s b^{(1)} + 2 \nu_s^{(1)} - \partial_0 c_s^{(1)}\right) . \end{aligned} \tag{E2c}$

      For $ X^{\mu} $ in Eq. (26), we have $ {\cal{X}}_{\mu \nu} $ in the form of

      $ \begin{aligned}[b] {\cal{X}}_{00}^S =& - 2 \left( 2 \phi^{(1)} + \frac{1}{a} \int \mathrm{d} \eta \left\{a \phi^{(1)} \right\} \partial_0 + \delta^{ij} \int \mathrm{d} \eta \left\{ \nu_i^{(1)} + \partial_i b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta'\left\{a \partial_i \phi^{(1)} \right\} \right\} \partial_j \right) \phi^{(1)} \\ & + 2 \delta^{ij} \left( \nu_j^{(1)} + \partial_j b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta \left\{a \partial_j \phi^{(1)} \right\} \right) \left(\partial_i b^{(1)} + \nu_i^{(1)}\right), \end{aligned}\tag{E3a} $

      $ \begin{aligned}[b] {\cal{X}}_{0 j}^S = & \Biggr( \delta^k_j\phi^{(1)} + \frac{\dot{a}}{a^2} \delta^k_j \int \mathrm{d} \eta\left\{a \phi^{(1)} \right\}+ \frac{1}{a} \delta^k_j \int \mathrm{d} \eta \left\{a \phi^{(1)} \right\} \partial_0 + \delta^{il}\delta^k_j \int \mathrm{d} \eta \Biggr\{ \nu_i^{(1)} + \partial_i b^{(1)} \\ & + \frac{1}{a}\int \mathrm{d} \eta' \left\{a \partial_i \phi^{(1)} \right\} \Biggr\} \partial_l + \delta^{lk} \left(\partial_j \int \mathrm{d} \eta \left\{ \nu_l^{(1)} + \partial_l b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta' \left\{a \partial_l \phi^{(1)} \right\} \right\}\right) \Biggr) \left(\partial_k b^{(1)} + \nu_k^{(1)}\right) \\ &+ \frac{2}{a} \phi^{(1)} \int \mathrm{d} \eta \left\{a \partial_j \phi^{(1)} \right\} + \delta^{ik} \left( \nu_i^{(1)} + \partial_i b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta \{a \partial_i \phi^{(1)} \} \right) \Biggr( - 4 \psi^{(1)} \delta_{jk} + 4 \partial_j \partial_k {e}^{(1)} \\ &+ 2 \partial_j c_k^{(1)} + 2 \partial_k c_j^{(1)} + 2 h_{j k}^{(1)} - \frac{2 \dot{a}}{a^2} \delta_{kj} \int \mathrm{d} \eta \left\{a \phi^{(1)} \right\} \\ &+ \left(\delta_{j}^s \partial_k + \delta_{k}^s \partial_j\right) \int \mathrm{d} \eta \left\{ \nu_s^{(1)} + \partial_s b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta' \{a \partial_s \phi^{(1)} \} \right\} \Biggr), \end{aligned}\tag{E3b} $

      $ \begin{aligned}[b] {\cal{X}}_{kl}^S =& \Biggr( \delta^s_k \delta^t_l \Biggr( \frac{2 \dot{a}}{a^2} \int \mathrm{d} \eta \{a \phi^{(1)} \}+ \frac{1}{a} \int \mathrm{d} \eta \{a \phi^{(1)} \} \partial_0 + \delta^{ij} \int \mathrm{d} \eta \Biggr\{ \nu_i^{(1)} + \partial_i b^{(1)} \frac{1}{a} \int \mathrm{d} \eta' \{a \partial_i \phi^{(1)} \} \Biggr\} \partial_j \Biggr) \\ &+ \left( \left(\delta^s_k \delta^{ti} \partial_l + \delta^t_l \delta^{s i} \partial_k\right) \int \mathrm{d} \eta \left\{ \nu_i^{(1)} + \partial_i b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta' \{a \partial_i \phi^{(1)} \} \right\} \right) \Bigg) \Biggr( - 4 \psi^{(1)} \delta_{st} + 4 \partial_s \partial_t {e}^{(1)} + 2 \partial_s c_t^{(1)} + 2 \partial_t c_s^{(1)} + 2 h_{s t}^{(1)} \\ &- \left(\delta_{s}^n \partial_t + \delta_{t}^n \partial_s\right) \int \mathrm{d} \eta \left\{ \nu_n^{(1)} + \partial_n b^{(1)} + \frac{1}{a} \int \mathrm{d} \eta' \{a \partial_n \phi^{(1)} \} \right\} - \frac{2 \dot{a}}{a^2} \delta_{st} \int \mathrm{d} \eta \{a \phi^{(1)} \} \Biggr) \\ &+ \frac{1}{a} \left(\partial_s b^{(1)} + \nu_s^{(1)}\right) \int \mathrm{d} \eta \Biggr\{a (\delta^s_k \partial_l + \delta^s_l \partial_k) \phi^{(1)} \Biggr\} . \end{aligned} \tag{E3c}$

      (130c)

    APPENDIX F: TRANSVERSE AND TRACELESS OPERATOR IN MOMENTUM SPACE
    • The transverse and traceless operator in the configuration space is expressed as

      $ \begin{aligned}[b] \Lambda^{lm}_{ij} = &{\cal{T}}^l_i {\cal{T}}^m_j - \frac{1}{2} {\cal{T}}_{ij} {\cal{T}}^{lm} \\ = &\left(\delta^l_i - \partial^l \Delta^{- 1} \partial_i\right) \left(\delta^m_j - \partial^m \Delta^{- 1} \partial_j\right)\\ &- \frac{1}{2} \left(\delta_{ij} - \partial_i \Delta^{- 1} \partial_j\right) \left(\delta^{lm} - \partial^l \Delta^{- 1} \partial^m\right) . \end{aligned} \tag{F1}$

      In momentum space, we could simply let $ \partial_l \rightarrow ik_l $. Eq. (F1) can be rewritten as

      $ \begin{aligned}[b] \Lambda^{lm}_{ij} =& \left( \delta^l_i - \frac{k^l k_i}{\left| {\boldsymbol{k}} \right|^2} \right) \left( \delta^m_j - \frac{k^m k_j}{\left| {\boldsymbol{k}} \right|^2} \right) - \frac{1}{2} \left( \delta_{ij} - \frac{k_i k_j}{\left| {\boldsymbol{k}} \right|^2} \right) \left( \delta^{lm} - \frac{k^l k^m}{\left| {\boldsymbol{k}} \right|^2} \right) \\ = &\left(\delta^l_i - n^l n_i\right) \left(\delta^m_j - n^m n_j\right)- \frac{1}{2} \left(\delta_{ij} - n_i n_j\right) \left(\delta^{lm} - n^l n^m\right), \end{aligned} \tag{F2}$

      where we have let $ n^l \equiv {k^l}/{\left| {\boldsymbol{k}} \right|} $ for simplicity.

      The three-dimensional momentum space can be spanned by a set of normalized orthogonal bases {$ {\epsilon}_i $, $ \bar{\epsilon}_i $, $ n_i $} that satisfy

      $ {\epsilon}_i {\epsilon}^i = \bar{\epsilon}_i \bar{\epsilon}^i = n_i n^i = 1, \tag{F3a}$

      $ {\epsilon}_i \bar{\epsilon}^i = {\epsilon}_i n^i = \bar{\epsilon}_i n^i = 0. \tag{F3b} $

      The indices are raised and lowered via the Kronecker delta $ \delta^{i}_{j} $, which can be written as

      $ \delta^j_i = {\epsilon}_i {\epsilon}^j + \bar{\epsilon}_i \bar{\epsilon}^j + n_i n^j . \tag{F4} $

      Substituting Eq. (F4) into Eq. (F2), we can obtain

      $ \begin{aligned}[b] \Lambda^{lm}_{ij} =& \left({\epsilon}_i {\epsilon}^l + \bar{\epsilon}_i \bar{\epsilon}^l\right) \left({\epsilon}_j {\epsilon}^m + \bar{\epsilon}_j \bar{\epsilon}^m\right) - \frac{1}{2} \left({\epsilon}_i {\epsilon}_j + \bar{\epsilon}_i \bar{\epsilon}_j\right) \left({\epsilon}^l {\epsilon}^m + \bar{\epsilon}^l \bar{\epsilon}^m\right) \\ =& {\epsilon}_i {\epsilon}_j {\epsilon}^l {\epsilon}^m + {\epsilon}_i \bar{\epsilon}_j {\epsilon}^l \bar{\epsilon}^m + \bar{\epsilon}_i {\epsilon}_j \bar{\epsilon}^l {\epsilon}^m + \bar{\epsilon}_i \bar{\epsilon}_j \bar{\epsilon}^l \bar{\epsilon}^m \\ &- \frac{1}{2} \left({\epsilon}_i {\epsilon}_j {\epsilon}^l {\epsilon}^m + {\epsilon}_i {\epsilon}_j \bar{\epsilon}^l \bar{\epsilon}^m + \bar{\epsilon}_i \bar{\epsilon}_j {\epsilon}^l {\epsilon}^m + \bar{\epsilon}_i \bar{\epsilon}_j \bar{\epsilon}^l \bar{\epsilon}^m\right) \\ =& \frac{1}{2} ({\epsilon}_i \bar{\epsilon}_j {\epsilon}^l \bar{\epsilon}^m + \bar{\epsilon}_i {\epsilon}_j {\epsilon}^l \bar{\epsilon}^m + \bar{\epsilon}_i {\epsilon}_j \bar{\epsilon}^l {\epsilon}^m + {\epsilon}_i \bar{\epsilon}_j \bar{\epsilon}^l {\epsilon}^m) \\ &+ \frac{1}{2} ({\epsilon}_i {\epsilon}_j {\epsilon}^l {\epsilon}^m - {\epsilon}_i {\epsilon}_j \bar{\epsilon}^l \bar{\epsilon}^m - \bar{\epsilon}_i \bar{\epsilon}_j {\epsilon}^l {\epsilon}^m - \bar{\epsilon}_i \bar{\epsilon}_j \bar{\epsilon}^l \bar{\epsilon}^m) \\ =& \frac{1}{\sqrt{2}} ({\epsilon}_i \bar{\epsilon}_j + \bar{\epsilon}_i {\epsilon}_j) \frac{1}{\sqrt{2}} ({\epsilon}^l \bar{\epsilon}^m + \bar{\epsilon}^l {\epsilon}^m) \\ &+ \frac{1}{\sqrt{2}} \left({\epsilon}_i {\epsilon}_j - \bar{\epsilon}_i \bar{\epsilon}_j\right) \frac{1}{\sqrt{2}} \left({\epsilon}^l {\epsilon}^m - \bar{\epsilon}^l \bar{\epsilon}^m\right) =: {\epsilon}^\times_{ij} {\epsilon}^{\times, lm} + {\epsilon}^+_{ij} {\epsilon}^{+, lm}, \end{aligned}\tag{F5} $

      where we define

      $ {\epsilon}^\times_{ij} \equiv \frac{1}{\sqrt{2}} ({\epsilon}_i \bar{\epsilon}_j + \bar{\epsilon}_i {\epsilon}_j), \tag{F6a} $

      $ {\epsilon}^+_{ij} \equiv \frac{1}{\sqrt{2}} ({\epsilon}_i {\epsilon}_j - \bar{\epsilon}_i \bar{\epsilon}_j). \tag{F6b} $

      Here, $ {\epsilon}_{ij}^+ $ and $ {\epsilon}^{\times}_{ij} $ are called the transverse and traceless polarization tensors. They are widely used to study the gravitational waves [61]. We note that the third equation in Eq. (F5) is established when $ \Lambda^{lm}_{ij} $ acts on the rank-two symmetric tensors.

Reference (62)

目录

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return