ABSTRACT (original article): We study second order gravitational perturbations on plane wave spacetimes from both the metric and curvature perturbation points of view. For the former, we explicitly use the isometries of the background to introduce tensor oscillator harmonics, which render Einstein’s equations algebraic around symmetric plane waves. For the latter, we formulate the first and second order Teukolsky equations in a Geroch-Held-Penrose covariant way. Both approaches are useful in their own right, and together with our discussion on gauge freedom, they provide a foundation for the study of higher-order gravitational dynamics around plane wave spacetimes. Taking the perspective that these plane wave spacetimes arise from Penrose limits, we subsequently use these results to explore the nonlinear gravitational dynamics close to black hole lightrings. Specifically, we define and discuss quadratic quasinormal mode ratios, observe that they satisfy emergent selection rules, and make publicly available a code to compute them. CITATION (original article): Fransen, K., Pereñiguez, D. & Redondo-Yuste, J. Perturbations of plane waves and quadratic quasinormal modes on the lightring. J. High Energ. Phys. 2025, 148 (2025). https://doi.org/10.1007/JHEP12(2025)148​
We will refer to this work as [1] or “the paper”.
Notebook: This notebook covers only the approach which we have found to be most computationally efficient in the paper: the curvature perturbation point of view represented by a second order Teukolsky equation. It is the definition and the extraction of the quadratic quasinormal mode ratios from this second order Teukolsky equation which we present here. At each step we will give some context and explanation. However, for a complete discussion of the theoretical background we refer to the paper.

0. Introduction & Overview

Quadratic quasinormal modes (QQNMs) represent the leading non-linear corrections to quasinormal modes – the characteristic resonances governing black hole relaxation to equilibrium. Quasinormal modes are already of interest in current gravitational wave observations and also QQNMs are expected to be observable in the future. They will provide a controlled, observational window into perturbative non-linearities of general relativity around black holes. ​For asymptotically flat black holes, quasinormal modes are in close relation with the bound null geodesics (the “lightring”). This correspondence can be made precise in the eikonal or high (real) frequency limit. Moreover, a version of the eikonal limit and the corresponding quasinormal modes can be captured by the physics of (simpler) plane wave spacetimes. It is then a natural question how such a correspondence extends to QQNMs and as one aspect of this question, the paper and the notebook below study QQNMs purely from a plane wave spacetime point of view. ​Optimistically, future work will relate these plane wave QQNMs in a precise, quantifiable way to a particular regime of QQNMs of real black holes. Regardless, they provide a beautiful toy model to help better understand such observable non-linearities in general relativity. ​This notebook computes: 1. Background structure (metrics, tetrads, spin coefficients, Weyl scalars) 2. Linear perturbations (first-order solutions) 3. Nonlinear perturbations (second-order solutions) 4. Quadratic QNM ratios (quantifying nonlinear effects)

1. Background Spacetime, Spin coefficients and GHP operators

The geometry of plane wave spacetimes is relatively simple and well - understood. We will not dwell on it long here and recommend readers to run the code and move on to the next section, perhaps glossing over “Some observations about the background” in order to be reminded of some key properties of the plane wave background and our chosen frame.
In the paper, the content roughly corresponds to Section 2 (Plane waves) and Appendix B (GHP formalism).

1.1. (Pseudo-)Riemannian geometry utilities

This subsection takes its routines directly from “EDCRGTcode.m”, a set of (pseudo-)Riemannian geometry routines originating from the work of S. Bonanos (see also https://library.wolfram.com/infocenter/ID/4484/). We reproduce the ones we need here explicitly in order to avoid external dependencies for user convenience. The main routine “RGtensors” takes a metric and a coordinate system to compute standard, associated operators and tensors such as the covariant derivative and the Riemann tensor. We refrain from further comments on these methods and recommend simply to run them and continue.
Quit[];

1.1.1. EDCRGTcode


1.2. Background spacetime

Next, we define the background plane wave spacetime into the variable “met” using the coordinates in “coords”. We define the associated covariant derivatives and curvature tensors using the utilities from previous section, the aforementioned “RGtensors” function. Remark that instead of immediately specializing to a plane wave spacetime, we first define a more general pp-wave spacetime using a general function H[u,x,y]. Later, we will use the replacement rule “SymmetricPlaneWaveRepl” to restrict ourselves to the particular homogeneous plane waves satisfying Einstein’s equations, studied in most detail in the paper.
coords={u,v,x,y};​​met=DiagonalMatrix[{-H[u,x,y],0,-1,-1}]+{{0,1,0,0},{1,0,0,0},{0,0,0,0},{0,0,0,0}};​​RGtensors[met,coords];​​SymmetricPlaneWaveRepl=FuncRepRules[H[u,x,y],-
2
Ω
(
2
x
-
2
y
)];​​​​$Assumptions={Ω∈Reals,Ω>0};
gdd =
-H[u,x,y]
1
0
0
1
0
0
0
0
0
-1
0
0
0
0
-1
LineElement = 2d[u]d[v]-
2
d[x]
-
2
d[y]
-
2
d[u]
H[u,x,y]
gUU =
0
1
0
0
1
H[u,x,y]
0
0
0
0
-1
0
0
0
0
-1
gUU computed in 0. sec
Gamma computed in 0. sec
Riemann(dddd) computed in 0.016 sec
Riemann(Uddd) computed in 0.015 sec
Ricci computed in 0. sec
Weyl computed in 0. sec
Einstein computed in 0. sec
All tasks completed in 0.03125 seconds

1.3. Spin coefficients and GHP operators


1.4. Some observations about the background


2. Linear Perturbations

We move on to describe the linear perturbations around plane wave backgrounds. As noted, we restrict ourselves to the curvature perturbation approach in this notebook. As discussed in Section 3 (Perturbation theory) of the paper, the Hertz potentials, from which linear metric perturbations can be constructed by the action of second derivative operators, is governed by the massless, scalar wave equation on the plane wave background.
​
For the symmetric plane waves, obtained here by the replacement defined in “SymmetricPlaneWaveRepl”, the solutions to the scalar wave equation with the correct asymptotic behavior to match onto eikonal quasinormal modes can be given analytically. After defining these in “2.1 Mode solutions”, with our particular normalization conventions defined in (4.7) of [1], we introduce various replacement rules for Hermite polynomials that will be useful as these polynomials feature in the solutions. Finally, we show that the defined modes are indeed solutions to the scalar wave equation.

2.1. Mode solutions


2.2. Utilities related to Hermite polynomials

In order to simplify the non-linear source and related second order perturbation solutions later, it will be useful to have available various identities that are satisfied by the Hermite polynomials and their products. These relations include for instance the recurrence relations as well as the expansion of products of Hermite polynomials into a sum of (single) Hermite polynomials. Such relations are fairly standard and can for instance be found in the digital library of mathematical functions: https://dlmf.nist.gov/ and in I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products (1943). See in particular (D.3) in the paper for the product of two Hermite polynomials. The “SumRule” below extracts only one term from the resulting sum, as specified by “dummy”.
MultRule[arg1_,arg2_,dummy_:r]:=HermiteH[a_,arg1]:>
a-2dummy
arg1
arg2
dummy
2
arg1
arg2
-1
Binomial[a,2dummy]
Factorial[2dummy]
Factorial[dummy]
HermiteH[dummy,arg2];​​xH=xHermiteH[a_,xb_]:>
1
2b
(HermiteH[a+1,xb]+2aHermiteH[a-1,xb]),x
2
HermiteH[a_,xb_]
:>
HermiteH[a,xb]
2b
(HermiteH[a+1,xb]+2aHermiteH[a-1,xb]),​​
2
x
HermiteH[a_,xb_]:>
4(-1+a)aHermiteH[-2+a,bx]+(2+4a)HermiteH[a,bx]+HermiteH[2+a,bx]
4
2
b
,
2
x
2
HermiteH[a_,xb_]
:>
1
4
2
b
(4(-1+a)aHermiteH[-2+a,bx]HermiteH[a,bx]+2a
2
HermiteH[a,bx]
+HermiteH[1+a,bx](2aHermiteH[-1+a,bx]+HermiteH[1+a,bx])),​​yHermiteH[a_,yb_]:>
1
2b
(HermiteH[a+1,yb]+2aHermiteH[a-1,yb]),y
2
HermiteH[a_,yb_]
:>
HermiteH[a,yb]
2b
(HermiteH[a+1,yb]+2aHermiteH[a-1,yb]),​​
2
y
HermiteH[a_,yb_]:>
4(-1+a)aHermiteH[-2+a,by]+(2+4a)HermiteH[a,by]+HermiteH[2+a,by]
4
2
b
,
2
y
2
HermiteH[a_,yb_]
:>
1
4
2
b
(4(-1+a)aHermiteH[-2+a,by]HermiteH[a,by]+2a
2
HermiteH[a,by]
+HermiteH[1+a,by](2aHermiteH[-1+a,by]+HermiteH[1+a,by]));​​SumRule[arg_,dummy_:]:=HermiteH[a_,arg]HermiteH[b_,arg]:>Binomial[a,dummy]Binomial[b,dummy]
dummy
2
Factorial[dummy]HermiteH[a+b-2dummy,arg],​​
2
HermiteH[a_,arg]
:>
2
Binomial[a,dummy]
dummy
2
Factorial[dummy]HermiteH[2a-2dummy,arg];

2.3. Check that modes are indeed solutions

Let us conclude our discussion of linear perturbations by showing that the modes defined in 2.1 indeed solve the linear wave equation, and that this is the Teukolsky equation on the plane wave spacetime as given in terms of GHP covariant derivatives in (2.16).
​
In order to act with the GHP covariant derivatives as defined previously, we first define the Hertz potentials with their GHP- weights
PsiH={ψH[u,v,x,y],-4,0};​​PsiHc={ψHc[u,v,x,y],0,-4};​​PsiHB={ψHB[u,v,x,y],-4,0};​​PsiHcB={ψHcB[u,v,x,y],0,-4};
Now we define the wave operator, and check its expression in coordinates.
box0[f_]:=2(þp[þ[PsiH]]-ðp[ð[PsiH]])[[1]]//.SymmetricPlaneWaveRepl//.FuncRepRules[ψH[u,v,x,y],f];​​box[f_]:=-D[f,{y,2}]-D[f,{x,2}]-(
2
x
-
2
y
)
2
Ω
D[f,{v,2}]+2D[f,u,v];​​box0[f[u,v,x,y]]==box[f[u,v,x,y]]//Simplify
True
Acting with this wave operator, we check that the general mode solutions are eigenfunctions for this operators with eigenvalue “eigvalue[pu,pv,nx,ny]”, that the “onShell[pv,nx,ny]” replacement rule restricts this to zero eigenvalues, and, hence, the “onShell$mode[pv,nx,ny]” modes are solutions to the free, massless scalar wave equation.
box[mode[pu,pv,nx,ny]]==eigvalue[pu,pv,nx,ny]*mode[pu,pv,nx,ny]//.xH//Simplify​​eigvalue[pu,pv,nx,ny]==0/.onShell[pv,nx,ny]//FullSimplify​​(box[onShell$mode[pv,nx,ny]]==0//.xH)//FullSimplify
True
True
True

3. Computation of Second Order Solutions

For the second order equations, our starting point is equation (3.40) in the paper together with the source terms (3.41) and (3.42), depending on the linear solutions. The derivations of these equations is one of the main technical results of the paper. However, our goal here is to solve these equations. Therefore, we simply start by defining the source terms as given in (3.41) and (3.42).
Subsequently, there are two steps in solving the second order equations
1
.
Re-expand the second-order source as a sum of first order modes (not on-shell) using the properties of the Hermite polynomials
2
.
Make use of the fact that the elementary modes are eigenfunctions of the linear wave operator to reduce the problem to an algebraic one for each term in this sum.
Step 1. will be the main objective of this section and is done in particular in “3.2. Routines to compute source term projections” below.
For more discussion, see also Section 4.2.2. in the paper as well as the illustrative toy model presented in Section 5.1 in the paper.
After step 1., as is discussed in Section 5.2 of the paper, it is in fact a matter of simple algebra to define the QQNM ratios, see in particular (5.18). We thus immediately continue to define the QQNM ratios and conclude this part of the notebook with some checks. Leaving an exploration of the QQNM ratios and their visualizations to the final part of the notebook.

3.1. Define source terms


3.2. Routines to compute source term projections


4. Compute Ratios

We now use the routines we have set up to compute the QQNM ratios to explore these across the parameter space of wave numbers. We will follow roughly Section 5.4 of the paper.

3.1. Coupling of Fundamental Modes

We start with the coupling of two “fundamental modes” with each other, meaning modes with nx=0 and ny=0. As explained in the paper, there are various qualitatively different behaviors depending on the relation between the relative values of pv between these modes. Specifically, there are the finite excitation channels, where only a finite number of quadratic modes are excited and the infinite excitation channels where an infinite number of quadratic modes are excited. The infinite number of modes is always in the mode numbers nx. In both cases, only a finite number of different ny modes can be excited. We illustrate both finite and infinite excitation channels below.
◼
  • First pa > pb > 0:
    ​
    For the + combination, this is a finite excitation channel and it only excites a fundamental quadratic mode:
    • To illustrate that no other modes are excited, for instance:
    For the - combination, this is an infinite excitation channel, although only even modes can be excited on account of the symmetries of the parent modes. We compare the result to the analytic expressions given in (5.31) in the paper.
    ◼
  • Next, pa > 0 > pb, pa+pb>0
    • For the - combination, this is a finite excitation channel and it only excites a fundamental quadratic mode:
    However, for the + combination, this is an infinite excitation channel, although again only even modes can be excited on account of the symmetries of the parent modes. We compare the result to the analytic expressions given in (5.30) in the paper.
    ◼
  • Finally, for the last two scaling limits of (5.33) in the paper, taking pa and pb comparably large although not identical in absolute value, note that
    • In summary, let' s visualize these self - couplings of the fundamental modes along the lines of Figure 1 in the paper. Note that we first combine the results from above in order to make the visualization slightly faster:

    3.2. Self-coupling of overtones

    Similarly, we can consider the self-coupling of modes more generally, so not necessarily fundamental but, for definiteness, with pa = pb
    In this case, to avoid zero modes, we do restrict ourselves to "RatioP", leaving us always with finite excitation channels when p>0 and maximal quadratic target mode of 2n:
    To illustrate that higher overtone target modes are not excited consider:

    5. Conclusion

    In the paper, we have studied perturbative non-linearities on plane wave spacetimes. Plane wave geometries are not only interesting in their own right, they also appear as approximate geometries near null geodesics. For instance, near equatorial light - rings of Kerr black holes, a special class of Lorentzian symmetric (Cahen-Wallach) plane waves emerge, which capture an asymptotic branch of quasinormal modes. The non - linearities around such plane waves can thus be interpreted and analyzed in terms of quadratic quasinormal modes. We therefore developed and used a perturbation framework to define and compute the excitation of quadratic quasinormal modes on such plane waves, and measure their excitation in the second-order fluctuations of the Weyl scalar Ψ0 in terms of QQNM ratios. In this notebook, we have provided you with the code to reproduce these QQNM ratios.

    CITE THIS NOTEBOOK

    Perturbations of plane waves and quadratic quasinormal modes on the lightring ​
    by Kwinten Fransen, David Pereñiguez, Jaime Redondo-Yuste
    Wolfram Community, STAFF PICKS, February 13, 2026
    ​https://community.wolfram.com/groups/-/m/t/3639726