Origin of Quasi-Periodic Oscillations and Accretion Process in X-Ray Binaries around Quantum Lee-Wick Black Hole

In recent years, there has been considerable interest in the perturbative approach to quantum gravity, particularly due to the new understanding of the challenges posed by higher derivatives and related ghost fields [1]. Several models and processes have been proposed to reconcile quantum gravity's unitarity with renormalizability. Interestingly, models of higher-derivative gravity with complex poles in the propagator, known as Lee-Wick gravity, have been investigated [2, 3]. The theory becomes super-renormalizable when higher derivatives are added to the action, improving the propagator's behavior in the ultraviolet limit [4]. As initially suggested by Lee and Wick [5, 6], the ghost-like poles that occur in complex conjugate pairs contribute to a unitary S-matrix [7−9].

Recent advances in understanding the problems of higher derivatives and the corresponding ghost fields have led to renewed interest in the perturbative approach to quantum gravity. To reconcile renormalizability and unitarity, two essential components of quantum gravity, researchers have put forth several models and processes. One prominent method is creating models of higher-derivative gravity, also known as Lee-Wick gravity [2, 3], which incorporate complex poles in the propagator. In other words, the propagator behaves better at high energies (also known as the ultraviolet regime) when higher derivatives are added to the action. A super-renormalizable theory [4] is the result of this enhancement. Moreover, ghost-like poles of complex conjugate pairs contribute to maintaining a unitary S matrix, guaranteeing that probabilities remain constant over time [2]. According to the relevant literature [7−9], this idea is reminiscent of the initial suggestions by Lee and Wick in their studies on quantum field theory [5, 6].

One of the most prominent features of complex poles at the classical level is the oscillating behavior seen in linearized solutions to the equations of motion [3, 10, 11]. These oscillations have been investigated both theoretically and experimentally and contribute to the gravitational force, especially in the low-energy regime [12−16]. However, higher derivatives complicate the analysis and make obtaining precise solutions to the entire set of field equations more challenging. To address this complexity, scientists have developed an empirical representation of Lee-Wick BHs that considers these oscillatory solutions [17]. This model aims to better understand how the oscillations affect BH behavior in the context of Lee-Wick gravity.

Recent research has investigated several important features of Lee-Wick BHs associated with the smeared source. The regularity of curvature and the thermodynamic properties of BHs [17−20, 22], the characteristics and conditions of the horizons for various sizes of BHs [20], the BHs impacted by noncommutative geometry [23, 24], gravitational light deflection [25−27], and the precession of orbits [28] are just a few of the topics covered in these investigations. The Newman-Janis algorithm [29] has also been used for a spinning extension of these BHs. An intriguing discovery, as stated in [20], is that the ratio $q = b/a$, where b and a reflect particular characteristics associated with the BH structure, affects the oscillation pattern of solutions. This ratio not only determines acceptable limits for the placement of the event horizon but also influences the possible number of horizons. In the end, the ratio q and the mass M of the effective source define the true number of BH horizons. A deeper understanding of the characteristics and behavior of Lee-Wick BHs under various physical circumstances may arise from this knowledge. Batic et al. [21] used the spectrum approach to investigate scalar, electromagnetic, and gravitational perturbations in their thorough analysis of quasinormal modes associated with Lee-Wick BHs. Our research takes into account all representative cases, both extremal and non-extremal, and shows that such BHs can exhibit a rich structure of horizons. Specifically, they demonstrated that for extreme and near-extremal configurations, entirely imaginary quasinormal modes appear, indicating a quick return to equilibrium without oscillation.

The sixth-derivative gravity field equations cannot be solved exactly by the Lee-Wick BHs suggested in [17, 20], but they might capture some of their key characteristics, namely curvature regularity and the existence of multiple or extreme horizons. The field equations derived from the extended Einstein-Hilbert action, containing all possible higher-derivative terms with up to six metric derivatives, have recently undergone a thorough investigation by researchers. They discovered that a bounded Kretschmann scalar $R_{\mu\nu\alpha\beta}^2$ [30] is present in all vacuum solutions that allow Frobenius expansion. Its equivalent to the Schwarzschild solution in general relativity appears to be a regular metric [30]. Current research and knowledge are limited to the weak field regime or deal with local features of the solutions [30, 31]. Due to this constraint, it is currently impossible to say with certainty if the horizon structures described in [20] match those of the precise solutions. However, it is noteworthy that the exact solutions of sixth-derivative gravity [30] contain both simple BHs and extreme horizons, confirming the requirement that a regular, asymptotically flat geometry must have an even number of horizons [31].

QPOs are rapid fluctuations in X-ray emissions from various astrophysical sources, such as Active Galactic Nuclei (AGN), high-mass X-ray binaries (HMXBs), low-mass X-ray binaries (LMXBs), and ultraluminous X-ray systems [32−34]. In LMXBs and HMXBs, BHs at their centers constantly accrete matter from a partner star [35, 37]. These oscillations are considered to result from the intricate interactions of matter spiraling inward into the BH. QPOs form as a result of these interactions, with oscillation frequencies ranging from millihertz to several hundred hertz [36]. The QPOs in AGNs can be discovered using measurements throughout the electromagnetic spectrum, including radio waves and gamma rays [38, 39]. However, it is crucial to note that this phenomenon is more commonly found in systems with stellar-mass BHs, as evidenced by X-ray spectrum data analysis.

The study of BHs has transcended the classical scope of general relativity to include aspects of quantum gravity, dark matter phenomenology, and horizon-scale studies. Regular BH models, which do not admit curvature singularities, have garnered enthusiasm for reconciling quantum spacetime geometry. Several works have shown how non-minimal couplings of scalars and curvature can modify Hawking BH radiation and modify the thermodynamics of BHs, as proposed by Garfinkle [40]. The well-established idea that certain classes of primordial or regular BHs, arising in modified gravity or quantum gravity-inspired scenarios, can act as viable candidates for dark matter, applies to both time-radial-symmetric and non-symmetric spacetime structures [41−43]. The dynamics of BHs and their associated gravitational-wave signatures have also been studied, particularly those signatures of dark-energy cosmological evolution [44]. Such ideas have encouraged studies of BH interactions with large-scale structures. The many revolutionary results from the Event Horizon Telescope have provided new avenues for studying BHs and their shadows to test fundamental physics theories. Horizon-scale imaging has been used to rule out black-hole solutions in mimetic gravity [45], impose the most detailed constraints on beyond-GR effects around Sgr A, and investigate the possible presence of ultralight bosonic particles via shadow distortions and superradiant instabilities. Moreover, studies on rotating regular BHs have uncovered significant connections between quasinormal modes and shadow features [46]. Other works have investigated the shadows of BHs with scalar hair [47], their role as standard candles in cosmology [48], and the shadow resulting from the presence of extra dimensions [49] or possible overspinning states of M87 [50].

Low-frequency quasi-periodic oscillations (LFQPOs) discovered from different sources have a frequency range from a few millihertz to many tens of hertz. To understand these frequencies, it may be necessary to examine the physical structures formed by matter accreting around BHs. This will help reveal important characteristics like the mass and angular momentum of the BHs at the cores of these systems. Observations from sources such as XTE J1550-564 have identified three different forms of LFQPOs: types A, B, and C [51−53]. Type-A QPOs have been discovered in about ten distinct sources and are quite uncommon in BH binary systems, typically appearing during the transition from the hard state to the soft state. Type-B QPOs are mostly associated with the physical processes behind astrophysical jets and have been observed in a limited number of BH binary systems. X-ray spectrum studies are typically used to detect Type-C QPOs, which are the most common in these systems [54−56]. GRS 1915+105 is an excellent example of a Type-C QPO and has been the focus of numerous theoretical and numerical investigations to clarify the physical properties of BH and the phenomena around it [53, 57−60]. The dynamic changes in matter flow around the BH's event horizon are thought to cause these oscillations.

The orbital dynamics of test particles in the framework of non-rotating Lee-Wick BH is the main topic of this publication. An outline of the paper's structure can be seen in the following: We present a comprehensive review of Lee-Wick BH in Section II. The ISCOs and equatorial circular orbits are discussed in the next section. We use an effective potential framework to examine the stability of these orbits. The oscillating motion of test particles around the Lee-Wick BH is examined in Section IV, which also covers the phenomenon of periastron precession. In addition, we present the dynamic structure of the shock cone produced by BHL accretion around a Lee-Wick BH using numerical simulation in Section V. The observational constraints of accretion regimes are presented in Section VI. In the last section, we provide the concluding remarks.

The effective field equations can be expressed as [17]:

where the Newton constant is denoted by G, and the Einstein tensor is represented by $ G_{\mu\nu} $. The effective stress-energy tensor becomes

It would replicate the effects of higher derivatives [17]. For a sixth-derivative Lee-Wick gravity model that is determined by the action,

With the conditions $ S_2 \gt 0,\; -2\sqrt{S_2} \lt S_1 \lt 2\sqrt{S_2} $ and the respective source leads to [17−19].

Here, the total source mass is denoted by M. With inverse length dimensions, the parameters a and b are connected to the couplings at work through

The real and imaginary parts of the Lee-Wick mass, $ m = a + i b $, are represented by a and b, respectively. Note that $ a, b \neq 0 $ and $ a, b \in \mathbb{R} $ are guaranteed by the domain of the parameters $ S_1 $ and $ S_2 $ in (3). By adding an effective equation of state and ensuring that the continuity equation $ \nabla_\mu \tilde{T}^\mu{}_\nu = 0 $ is satisfied [17, 20, 22], the remaining components of $ \tilde{T} $ in Eq. (2) can be identified. This method is comparable to that used for black holes (BHs) affected by noncommutative geometry [23, 24] and other effective metrics [19]. The line element describing the geometry of a spherically symmetric Lee-Wick BH is given by [21]:

where the lapse function $ f(r) $ reads:

Here, M is the mass of the gravitational source. The parameters $ S_1 $ and $ S_2 $ are related to the so-called Lee-Wick "mass" by $ \tilde{m} = S_1 + i S_2 $ with $ S_1, S_2 \gt 0 $. The rationale underlying the choice of signs for the parameters $ S_1 $ and $ S_2 $ can be found in [21]. It is noteworthy that, in the present work, we consider parameter values both within the admissible region, specified by $ S_2 \gt 0 $ and $ -2\sqrt{S_2} \lt S_1 \lt 2\sqrt{S_2} $, and outside this constraint set, in order to investigate the influence of these two parameters on our analysis. Horizons exist at those values of r where $ f(r) $ not only vanishes but also changes sign when crossing these points. Furthermore, the case of a Schwarzschild black hole is recovered in the asymptotic regime $ r \gg 1/S_1 $. The behavior of the lapse function $ f(r) $ as a function of r for different values of the parameters $ S_1 $ and $ S_2 $ is shown in Fig. II. The first column is plotted for $ S_2 = 1 $, and varying the parameter $ S_1 $. The second and third columns are plotted for fixed $ S_1=1 $ and $ S_1=2 $, respectively, while varying $ S_2 $. We see in all cases that the function $ f(r) $ crosses the axis $ f(r)=0 $. This means that for these ranges of parameters, the black hole solution exists. For the parameter $ S_1=1 $, the lapse function changes as the parameter $ S_2 $ changes with increasing radius r. However, for the parameter $ S_1=2 $, the lapse function does not change as the parameter $ S_2 $ changes with increasing radius r.

Figure 2.  (color online) Plots of the effective potential of particles moving around a non-rotating Lee-Wick BH.

Figure 1.  (color online) Below is the revised text: Behavior of the lapse function of the spacetime of a non-rotating Lee-Wick black hole.

The motion of a neutral particle can be described by the Hamiltonian given by

where μ is the mass of the particle, $ p^{\alpha}=\mu u^{\alpha} $ represents the four-momentum, $ u^{\alpha}= dx^{\alpha}/d\tau $ denotes the four-velocity, and τ is the proper time of the test particle. The Hamilton equations of motion can be written as

where $ \zeta=\tau/\mu $ denotes the affine parameter. Owing to the symmetries inherent in the black hole (BH) spacetime geometry, there exist two conserved quantities along geodesic motion, namely the specific energy $ {\cal{E}} $ and the specific angular momentum $ {\cal{L}} $, which are given by $ \dfrac{p_{t}}{m} = g_{tt} u^t + g_{t\phi} u^\phi = -{\cal{E}},$ where $ {\cal{E}}=E/m $, $ {\cal{L}}=L/m $ are the specific energy and specific angular momentum. The time $ u^t $, azimuthal $ u^\phi $, and radial $ u^r $ components of the four-velocity $ u^\alpha $ satisfy the following equations of motion:

where $ \epsilon = 1 $ for timelike particles and $ \epsilon = 0 $ for null-like particles. The dot denotes the derivative with respect to proper time τ. In our work, we focus only on timelike particles. The Hamiltonian by Eq. (III), for a non-rotating Lee-Wick BH, can be written in the form

where the effective potential $ V_{eff}(r, \theta) $ takes the form:

The effective potential $ V_{eff}(r, \theta) $ is crucial for understanding the motion of test particles, as it allows one to describe the trajectories of the particles without directly solving the equations of motion. The maximum and minimum values of $ V_{eff} $ correspond to unstable and stable circular orbits, respectively. Fig. 2 illustrates how $ V_{eff} $ varies with r. Stable circular orbits are associated with the minima, whereas unstable orbits occur at the maxima of the potential. The circular orbits for the equatorial plane are found by the following conditions:

To solve Eq. (III A), we obtain the circular orbits around a non-rotating Lee-Wick BH.

and the corresponding energy function takes the form:

Fig. 3 shows the angular momentum $ {\cal{L}} $ of equatorial circular orbits around a non-rotating Lee-Wick BH. The first column illustrates the influence of the parameter $ S_1 $ on the radial profiles of angular momentum. We observe that there is not much difference in the angular momentum while varying the parameter $ S_1 $. The second and third columns are plotted for $ S_1=1 $ and $ S_1=2 $, respectively, varying the parameter $ S_2 $. For $ S_1 = 1 $, the angular momentum decreases with increasing values of the parameter $ S_2 $ until $ r=6 $, then increases in the range $ (6,8) $. However, when $ S_1=2 $, the parameter $ S_2 $ does not change dramatically and almost exhibits the same behavior as the $ S_1 = 1 $ case.

Figure 3.  (color online) Plots of angular momentum of particles moving around a non-rotating Lee-Wick BH.

Fig. 4 presents the energy $ {\cal{E}} $ of equatorial circular orbits around a non-rotating Lee-Wick BH. The first panel shows how parameter $ S_1 $ affects the radial profiles of energy, while the second and third panels illustrate the impact of parameter $ S_2 $. As the parameter $ S_1 $ increases, the energy of the particles decreases. An increase in the parameter $ S_2 $ results in a decrease in particle energy for a given $ S_1 $.

Figure 4.  (color online) Plots of the energy of particles moving around a non-rotating Lee-Wick BH.

The locations of stable and unstable circular orbits correspond to the minimum and maximum of the effective potential, respectively. In Newtonian theory, the effective potential always has a minimum for any given angular momentum, and there are no ISCOs with a minimum radius. However, when the form of the effective potential depends on the angular momentum of the particle and other parameters, the situation changes. In general relativity, for particles orbiting near a Schwarzschild black hole, the effective potential exhibits two extrema for any angular momentum value. Only for a specific value of the angular momentum do these two points coincide, marking the ISCO, which is located at $ r = 3r_g $, where $ r_g $ is the Schwarzschild radius. The ISCO can be determined by applying certain conditions:

One can obtain the radii of ISCOs from the last relation around a non-rotating Lee-Wick BH. Fig. 5 illustrates the equatorial ISCOs around a non-rotating Lee-Wick BH. The first panel shows how the ISCOs change for different values of the parameter $ S_1 $, while the second panel depicts the influence of the parameter $ S_2 $. As the value of the parameter $ S_1 $ increases, the radii of ISCOs initially decrease and then become constant. Similarly, increasing the parameter $ S_2 $ leads to an initial decrease in the radii of the ISCOs, after which they become constant.

Figure 5.  Plots for ISCOs of particles moving around a non-rotating Lee-Wick black hole (BH).

To study the oscillatory motion of neutral particles, we perturb the equations of motion around stable circular orbits. When a test particle is slightly displaced from its equilibrium position associated with a stable circular orbit in the equatorial plane, it exhibits epicyclic motion, defined by linear harmonic oscillations.

The frequencies of harmonic oscillatory motion, as determined by the local observer, are given by $ \omega_{r}^{2} = -\dfrac{1}{2} \dfrac{\partial^{2} {\rm{V_eff}} (r, \theta)}{\partial r^{2}},$ The radial $ (\omega_{r}) $, latitudinal $ (\omega_{\theta}) $, and orbital/axial $ (\omega_{\phi}) $ frequencies of the neutral test particle for a non-rotating Lee-Wick BH take the form

where

The locally determined angular frequencies $ \omega_\beta $ are specified in Eqs. (17)-(19). In contrast, the angular frequencies measured by a static distant observer (Ω) are given by $\Omega_{i} = \omega_{i} \dfrac{{ \mathrm{d}} \tau}{{ \mathrm{d}} t},$

where ${ \mathrm{d}} \tau/{ \mathrm{d}} t$ is the redshift coefficient, given by

If the frequencies of small harmonic oscillations are measured in physical units by a distant observer, then the frequencies of neutral particles determined by these distant observers take the form

where $i\in\{r,\theta,\phi\}$; $\Omega_{r}$, $\Omega_{\theta}$, and $\Omega_{\phi}$ denote the dimensionless radial, latitudinal, and axial angular frequencies measured by a distant observer, given in the Appendix.

Figs. 6 and 7 present the radial profiles of the frequencies $\nu_{j}$ for small harmonic oscillations of neutral particles around a non-rotating Lee-Wick BH, as observed by a distant observer. Different cases of parameters are plotted. The radial profiles of the orbital frequency ($\Omega_{\phi}$) and the latitudinal frequency ($\Omega_{\theta}$) coincide. As the BH parameter $S_2$ increases, the profiles shift closer to the event horizon of the BH. In particular, an increase in the parameter $S_1$ also causes the profiles to move towards the BH horizon.

Figure 6.  (color online) Plots of the frequencies of particles moving around a non-rotating Lee-Wick BH with $\theta = \pi/2$ show that $\nu_{\theta} = \nu_{\phi}$. Consequently, the corresponding curves exhibit identical behavior.

Figure 7.  (color online) Plots of frequencies of particles moving around a non-rotating Lee-Wick BH with $\theta = \pi/2$, which implies $\nu_{\theta} = \nu_{\phi}$. Consequently, the corresponding curves exhibit identical behavior.

In this section, we examine the periastron frequency of a neutral test particle orbiting a non-rotating Lee-Wick BH, focusing on small perturbations around the equatorial plane at $ \pi/2 $. To calculate the periastron precession, we consider a slight perturbation of the particle from its stable position, which leads to oscillations around that stable point characterized by a radial frequency $ \Omega_{r} $. The periastron frequency, denoted $ \Omega_{P} $, is defined as the difference between the orbital frequency $ \Omega_{\phi} $ and the radial frequency $ \Omega_{r} $, as expressed in the following relation:

Unlike Newtonian gravity, where the radial and orbital frequencies are equal, general relativistic effects near a BH lead to the inequality $ \Omega_{\theta} \neq \Omega_{\phi} $. Fig. 8 shows the graphical behavior of the periastron frequency for particles orbiting a non-rotating Lee-Wick BH, presenting the effects of different BH parameters $ S_1 $ and $ S_2 $. The first column is plotted for different values of the parameter $ S_1 $. As the parameter $ S_1 $ increases, the frequency of the periastron increases. The second and third columns are plotted for different values of the BH parameter $ S_2 $. For $ S_1=1 $, different radial profiles can be observed for different values of the BH parameter $ S_2 $. However, for $ S_1=2 $, different profiles coincide.

Figure 8.  (color online) Plots for the periastron frequency of particles around a non-rotating Lee-Wick black hole (BH).

In this section, in addition to the theoretical developments, we investigate the dynamic structure of the shock cone produced by BHL accretion around a Lee-Wick BH, as well as the modifications induced by Lee-Wick gravity, by numerically solving the general relativistic hydrodynamics equations. To this end, we adapt the numerical code used in our previous work [61−67], configuring it with the appropriate initial and boundary conditions. The plasma distribution and the shock cone that form around the Lee-Wick BH are initialized using the parameter sets provided in Table 1. Each row of Table 1 lists the model name, the Lee-Wick parameters $S_1$ and $S_2$, the inner and outer horizons of the BH for that model, and the inner radial boundary of the computational domain (i.e., the radius of the innermost disk closest to the BH horizon). In the first block of Table 1, we model cases where the event horizon is close to $2M$ or slightly smaller. Consequently, the flow yields Schwarzschild-like outcomes, and below we present numerical results together with comparisons to the Schwarzschild case. In contrast, the second block considers cases with an event horizon radius $>2M$, as illustrated in Fig. 2, especially for $S_1=1$, where the altered behavior of the lapse function manifests in plasma morphology and in the formation of the cone and bow shock.

In the non-rotating Lee-Wick BH metric given by Eqs. II and 7, the lapse function $f(r)$ depends on the parameters $S_1$ and $S_2$. The parameter $S_1$ sets the exponential decay scale of the Lee-Wick correction, while $S_2$ controls its oscillatory behavior. As shown in Fig. 2, the zeros of the lapse, $f(r)=0$, determine the inner and outer horizons, which are listed in Table 1 for the corresponding pairs $(S_1,S_2)$. In the metric, one finds that for $r \gg 1/S_1$, the spacetime approaches the Schwarzschild limit, and the sensitivity of $f(r)$ to $S_2$ diminishes as $S_1$ increases (as seen in the lapse-function sweeps). In contrast, for $S_1=1$, there is a pronounced variation with $S_2$, whereas for $S_1=2$, almost no change is observed since the event horizons form at nearly the same radii as in the Schwarzschild case. Therefore, the choice of $(S_1,S_2)$ significantly modifies the shock cone morphology produced by BHL accretion. Oriented by the lapse behavior in Fig. 2 and the two initial data blocks listed in Table 1, we performed the numerical simulations accordingly.

It should be emphasized that the numerical models listed in Table 1 were deliberately selected from parameter sets lying both inside and outside the admissible region defined by $S_2>0$ and $-2\sqrt{S_2}<S_1 <2\sqrt{S_2}$. This choice allows us to investigate not only the weak- and strong-correction regimes of Lee-Wick gravity but also the different horizon structures that arise in the spacetime, including regular black-hole configurations, multiple-horizon cases, and horizonless/naked-singularity solutions. Therefore, the numerical analysis is designed to reveal how the accretion flow morphology, shock-cone structure, and QPO signatures change across both admissible and non-admissible sectors of the parameter space.

In the Block-1 models presented in Table 1, only the values of $ S_{1} $ and $ S_{2} $ that produce weak correction terms around the Lee-Wick BH are considered. These parameters lead to only small deviations from the classical Schwarzschild geometry. In this domain, the exponential suppression controlled by $ S_{1} $ dominates over the oscillatory contribution generated by $ S_{2} $, so the metric describing Lee-Wick gravity preserves the essential features of general relativity and remains nearly unchanged. Only higher-derivative terms introduce mild modifications in the vicinity of the BH horizon. The spacetime metric remains smooth and regular, and the position of the horizon exhibits only minimal shifts compared to the Schwarzschild case. Therefore, this region corresponds to the near-classical limit of the Lee-Wick model, providing a consistent setting for examining how subtle quantum-correction terms modify the spacetime geometry. In this section, we reveal how these small corrections influence the morphology of the shock cone and the characteristics of the QPOs obtained from the numerical simulations.

Understanding how the morphology of the shock cone formed in the downstream region of the BHL accretion flow changes with the Lee-Wick parameters $ S_{1} $ and $ S_{2} $ around a nonrotating BH helps reveal their influence on both the accretion structure and the resulting QPOs. For this purpose, Fig. 9 presents the dynamic structure of the shock cone obtained for the Block-1 parameter sets, together with the corresponding case of a Schwarzschild BH for direct comparison. The figure shows the rest-mass density distribution of the accreting matter on the equatorial plane. Both contour and vector plots are included for each model, allowing the behavior of the flow to be examined once the cone reaches a steady state. In all cases, the system clearly attains a quasi-steady configuration. The weak correction regime, where $ S_{1} $ and $ S_{2} $ introduce only mild deviations from the Schwarzschild spacetime, yields an overall shock cone morphology that is almost identical to the classical BHL configuration. The cone formed in the downstream region is smooth, symmetric, and tightly bound to the inner boundary of the computational domain at $ r = 2.3M $. This demonstrates that the plasma dynamics are primarily governed by the standard general relativistic potential rather than by the higher derivative corrections of Lee-Wick gravity. However, as $ S_{2} $ increases, the density contours along the cone begin to show slight asymmetry, and the opening angle of the cone becomes marginally smaller. The dependence of this angle on the parameters, shown in Fig. 11 and discussed later in detail, is consistent with the exponential oscillatory competition appearing in the lapse function $ f(r) $ given by Eq. 7. The numerically obtained shock morphologies in the weak Lee-Wick regime agree fully with the theoretical expectations derived from the effective potential profiles presented in Fig. 2. In particular, larger values of $ S_{1} $ produce strong exponential suppression in $ f(r) $, deepening the gravitational potential well and creating a more compact, high-density region near $ r \simeq 2M $. Meanwhile, increasing $ S_{2} $ induces oscillatory modulation in both $ f(r) $ and the curvature terms, slightly perturbing the equilibrium structure without destabilizing the flow trapped within the cone.

Figure 9.  (color online) In the Block-1 case given in Table 1, together with the Schwarzschild solution, we plot the rest-mass density in the region close to the BH, along with its contour map and vector (velocity) field on the equatorial plane, to display the morphology of the resulting shock cone.

Figure 11.  In the Block-1 configuration, we present the dependence of the shock cone opening angle on $ S_1 $ and $ S_2 $, shown alongside the Schwarzschild case. Consistent with theoretical expectations, at sufficiently large $ S_1 $, the opening angle converges to the Schwarzschild result.

In the weak-correction regime, the azimuthal variation of the rest-mass density inside the shock cone formed by the BHL accretion mechanism around the nonrotating Lee-Wick BH is illustrated in Fig. 10. The top panel shows this variation at $ r = 2.66M $, very close to the BH horizon, while the bottom panel presents the corresponding distribution near the ISCOs region. The rest-mass density profiles at both radial points are compared with the Schwarzschild case to demonstrate the dependence of the shock cone morphology on the Lee-Wick parameters $ S_1 $ and $ S_2 $ within the weak field regime. At $ r = 2.66M $, the cone is compact and well-developed. The slight deviations observed for different combinations of $ (S_1, S_2) $ arise from the modifications of gravity induced by the Lee-Wick corrections. As the parameter $ S_1 $, which controls the exponential suppression term in the lapse function $ f(r) $, increases, the resulting deeper effective potential well leads to a denser and more confined plasma region near the cone axis. In contrast, an increase in $ S_2 $, which regulates the oscillatory component of $ f(r) $, generates small azimuthal fluctuations and introduces a slight asymmetry compared to the Schwarzschild case. The insets inside the top panel clearly display the dependence of both the maximum density and the width of the strong shock region on the Lee-Wick parameters. For large values of $ S_1 $ ($ S_1 \geq 2 $), these variations converge towards the Schwarzschild solution, indicating that strong exponential damping suppresses the higher derivative corrections.

Figure 10.  (color online) It shows the dynamic structure of the shock cone formed around the non-rotating Lee-Wick BH and the parameter-dependent variation of the rest-mass density trapped inside the cone along the azimuthal direction at two different radial positions. In the upper panel, the variation at $ r = 2.66M $, the closest point to the BH, is presented for different values of $ S_1 $ and $ S_2 $ given in Block-1, in comparison with the Schwarzschild case. The small insets within the figure display the dependence of the strong shock region and the maximum density on the parameters. In contrast, the lower panel shows how the maximum density of matter trapped inside the cone varies with the parameters at $ r = 6.11M $, i.e., near the ISCOs, also compared with the Schwarzschild case.

In the bottom panel of Fig. 10, the azimuthal variation of the rest-mass density is shown at $ r = 6.11M $, corresponding to the region where the flow transitions from purely radial accretion to rotational motion near the ISCOs. Compared with the $ r = 2.66M $ case, the density amplitude decreases in all models due to the weaker gravitational confinement at larger radii. However, the dependence on $ S_1 $ and $ S_2 $ follows the same trend observed in the top panel. At $ r = 6.11M $, larger values of $ S_1 $ tend to make the azimuthal density distribution inside the cone smoother and more concentrated around the cone axis, whereas larger $ S_2 $ values induce fine-scale oscillatory modulations in the density profile. We emphasize that this statement refers to the local density structure at a fixed radius, not to a monotonic change in the global shock-cone opening angle. This region is particularly important since QPOs can emerge here, linking the spatial morphology of the cone to the temporal fluctuations in the accretion rate, as shown later in Fig. 12. Compared to the Schwarzschild case, both radial positions reveal that the Lee-Wick parameters $ S_1 $ and $ S_2 $ produce observable modifications in the shock cone structure while maintaining a symmetric and smooth overall geometry. The numerical results obtained here are in excellent agreement with the theoretical predictions derived earlier in the paper.

Figure 12.  (color online) Time evolution of the mass accretion rate at $ r = 2.3M $ around the non-rotating Lee-Wick BH for the Block-1 models. The left panel compares the Schwarzschild case with models at fixed $ S_1 = 6 $ and varying $ S_2 $, while the right panel shows models at fixed $ S_2 = 1 $ and varying $ S_1 $. The figure illustrates the strong instabilities of the plasma flow, the resulting shock-cone oscillations, and the dependence of the accretion-rate variability on the Lee-Wick parameters.

In addition to Fig. 10, Fig. 11 illustrates how the opening angle of the shock cone varies with the Lee-Wick parameters. As seen in Fig. 11, increasing $ S_1 $ causes the cone opening angle to gradually converge to the value obtained for the Schwarzschild solution. This confirms that the strong exponential suppression term minimizes deviations from general relativity, consistent with the theoretical discussion above. Meanwhile, variations in $ S_2 $ produce noticeable changes in cone width, reflecting the interaction between oscillatory and exponential terms in the lapse function. Fig. 11 shows that the global opening angle of the shock cone is not a monotonic function of $ S_1 $. The widest cone appears for the horizonless case $ S_1=1, S_2=1 $, where the weaker confinement broadens the downstream flow. For intermediate parameter values, the opening angle decreases, indicating stronger gravitational focusing and a more compact inflow geometry. For sufficiently large $ S_1 $, however, the opening angle converges toward the Schwarzschild value, consistent with the strong exponential suppression of the Lee-Wick corrections. This result is consistent with the improved compression and narrow cone geometry shown in Fig. 10. Altogether, these findings confirm that even mild quantum gravity corrections can leave subtle yet observable imprints on the accretion dynamics around the Lee-Wick BH.

When considering the case of $ S_1 = 1 $ and $ S_2 = 1 $ in Fig. 11, we observe analytically that no classical BH horizon is formed, as seen in Table 1. In this configuration, the spacetime corresponds to a naked singularity. Consequently, as the accreting matter falls toward the center, the gravitational pull acting on it at $ r = 2.3M $ is significantly weaker compared to the regular BH cases. Due to this reduced gravitational acceleration, the inflowing plasma becomes less focused and tends to spread more before reaching the downstream region. As a result, a shock cone is formed with a wider opening angle, lower density, and weaker compression. This broader and less confined cone structure is a direct manifestation of the weaker gravitational confinement around the naked singularity.

In the previous analysis, we demonstrated the formation of a shock cone around the nonrotating Lee-Wick black hole (BH) and showed that its dynamic structure varies with the model parameters. In this section, we calculate the mass accretion rate to investigate the possible emergence of quasi-periodic oscillations (QPOs). For this purpose, the accretion rate is measured at $ r = 2.3M $, corresponding to the inner boundary of the computational domain, where the inflowing matter leaves the simulation grid and falls into the BH. By evaluating the time-dependent behavior of $ dM/dt $, we determine how the accretion rate changes concerning the Lee-Wick parameters $ S_1 $ and $ S_2 $. Fig. 12 shows the temporal evolution of the mass accretion rate at $ r = 2.3M $ for the Block-1 models, which include different combinations of $ (S_1, S_2) $, along with the Schwarzschild reference case for direct comparison. All models exhibit noticeable temporal variations in the accretion rate, confirming that the accretion flow and the shock cone surrounding the BH are inherently dynamic rather than completely stationary. In the Schwarzschild case, the accretion rate shows small but distinct periodic oscillations, demonstrating that even without any quantum corrections, the downstream shock cone produces moderate and self-sustained instabilities. However, these oscillations remain relatively regular and have lower amplitudes. In contrast, Lee-Wick BH models show stronger and less regular variability in $ dM/dt $, with the amplitude and frequency of oscillations showing a mild dependence on the parameters $ S_1 $ and $ S_2 $.

As the parameter $ S_1 $ increases, the exponential suppression term in the lapse function modifies the effective potential, resulting in a denser inflow of matter and higher amplitude oscillations. On the other hand, an increase in $ S_2 $ enhances the oscillatory nature of the metric, producing more complex and quasi-periodic modulations that differ from the smoother Schwarzschild pattern. This parameter-dependent behavior indicates that the interplay between the exponential and oscillatory corrections in the Lee-Wick metric destabilizes the post-shock flow, leading to turbulence-like variability and the emergence of QPO-like features. Although the Schwarzschild solution already exhibits intrinsic oscillations associated with the breathing of the shock cone, the Lee-Wick configurations amplify these oscillations and modify their characteristic frequencies. The degree and nature of these changes are discussed in detail in Section V A 2. The time-dependent behavior of the accretion rate, particularly for large values of $ S_2 $, confirms the formation of standing shock wave structures and strong global oscillation modes within the downstream region. These enhanced oscillatory modes could represent potential observational signatures of Lee-Wick gravity. Although these results are obtained under weak-field initial conditions, the differences observed when compared to the Schwarzschild solution suggest that such variability patterns may offer a possible avenue for testing Lee-Wick gravitational effects in future observational or numerical studies.

To test the physical detectability of the Lee-Wick BH, it is essential to investigate the QPO modes that are trapped and excited within the shock cone mechanism forming around the BH. A comparative analysis between the Lee-Wick and Schwarzschild spacetimes allows us to identify both their similarities and differences, revealing how observational signals could carry the imprints of Lee-Wick gravity. Based on this motivation and considering that the oscillatory patterns shown in the two panels of Fig. 12 arise from the variability of the mass accretion rate, we performed a PSD analysis of the accretion rate. Fig. 13 and Fig. 14 present the PSD computed using the mass-accretion rate at $ r = 2.3M $, for various combinations of parameters $ (S_{1}, S_{2}) $ of the Schwarzschild and Lee-Wick models. To express the resulting frequencies in units of Hertz, we set the BH mass at $ M = 10 M_{\odot} $. After confirming that the QPO modes that arise in the PSD are global in nature, rather than the product of local instabilities, the analysis was carried out using only the accretion rate at $ r = 2.3M $. Numerical experiments performed with the accretion rate at $ r = 6.11M $ produced identical QPO frequencies, confirming the robustness of the signal. In Fig. 13, the PSD derived from the strongly oscillatory accretion rate at $ r = 2.3M $ displays the QPO peaks around the Schwarzschild and Lee-Wick BHs for constant $ S_{1} $ and varying $ S_{2} $. The upper left panel corresponds to the Schwarzschild case, showing a single broad LFQPO peak, which represents the typical post-shock oscillation observed in classical BHL accretion. When the Lee-Wick correction term modifies the spacetime geometry, additional sharp peaks appear in the PSD. The frequencies of these peaks increase systematically with larger $ S_{2} $, causing the QPOs to shift from the LFQPO to the high-frequency QPO (HFQPO) regime. This result, consistent with theoretical expectations, confirms that the strong oscillatory term in the Lee-Wick metric amplifies the curvature-driven instabilities within the shock cone.

Figure 13.  (color online) The Power Spectral Density analysis is performed at $ r = 2.3M $ using the mass accretion rate in Fig. 12 to reveal the quasi-periodic behavior around both the Schwarzschild and Lee-Wick BHs for different values of $ S_1 $ and $ S_2 $ as given in Block-1. In the top-left panel, the PSD analysis and the corresponding QPO frequencies are shown only for the Schwarzschild case, while in the other panels, both Schwarzschild and Lee-Wick BHs are plotted together. However, the peak frequencies indicated in the PSD correspond only to those that emerge in the Lee-Wick BH cases. As clearly seen from the PSD analyses, increasing values of $ S_2 $ cause the QPO frequencies to shift from LFQPO to HFQPO regimes.

Figure 14.  (color online) The same analysis as in Fig. 13 is performed, but this time the PSD results are shown for varying $ S_1 $ while keeping $ S_2 $ fixed. Similar quasi-periodic behavior is also observed in this case.

For low values of $ S_{2} $ ($ S_{2}=1 $), the observed peaks exhibit an approximate 3:2 ratio, likely arising from nonlinear resonances between the radial and azimuthal oscillation modes. As $ S_{2} $ increases, these ratios evolve toward 2:1 or 3:1, revealing higher harmonics and complex mode coupling indicative of stronger Lee-Wick effects. Moreover, similar mode ratios are observed in different combinations of $ (S_{1},S_{2}) $, consistent with observational data. The rich set of peaks obtained in each model arises from the superposition of fundamental oscillation modes, leading to new resonant frequencies. Consequently, when $ S_{1} $ is kept constant, increasing $ S_{2} $ increases the probability that the observable QPO frequencies correspond to HFQPOs.

In Fig. 13, it was shown that the increase in the parameter $ S_2 $ in the PSD analysis causes the characteristic frequencies to evolve toward the HFQPO regime. Consequently, the mode ratios emerging from the simulations are observed to be most prominently $ 3:2 $, with additional branches $ 2:1 $ and $ 4:3 $. In particular, the numerically detected $ 3:2 $ HFQPO pairs correspond well to those observed in stellar-mass BH systems, such as $ GRO J1655-40 $ ($ 300/450 $), $ XTE J1550-564 $ ($ 184/276 $), and $ GRS 1915+105 $ ($ 41/67 $) [68−71]. These observational findings confirm the presence of stable $ 3:2 $ frequency ratios in real astrophysical systems. However, the appearance of $ 2:1 $ ratios corresponds to transient states typically observed during state transitions in X-ray binaries. Consequently, the numerically calculated QPO ratios around the Lee-Wick BH, especially for parameter sets $ S_1 = 6 $, $ S_2 = 1 $ (clear $ 3:2 $) and $ S_1 = 6 $, $ S_2 = 6 $ ($ \approx 3:2$ with a nearby $ 2:1 $ branch), are in strong agreement with the observed HFQPO ratios in BH binaries. The presence of additional $ 4:3 $ branches reflects intermediate resonant behavior, which has also been discussed in relativistic resonance interpretations of QPO phenomena.

In Fig. 14, the PSD analysis is repeated, keeping $ S_{2} $ fixed and varying $ S_{1} $. In this case, the quasi-periodic nature of the signal persists across all models, but the dominant peaks gradually shift toward lower frequencies, indicating a transition toward LFQPOs. This behavior is attributed to the strong exponential damping term, which deepens the effective potential and lengthens the characteristic oscillation period of the accretion flow. Due to nonlinear coupling among oscillation modes, the PSDs exhibit a rich structure, making the 3:2 and 2:1 frequency ratios highly likely to be observed around the Lee-Wick BH. As seen in Fig. 14, when $ S_{2} $ is fixed, increasing $ S_{1} $ not only modifies the frequency positions but also enhances the amplitude of the PSD peaks. Thus, the detectability of such peaks in observational data becomes more probable.

For the models presented in Fig. 14, the mode ratios obtained from numerical calculations are highly consistent with the observed HFQPO frequency pairs in X-ray binary systems. For the case of $ S_1 = 1 $ and $ S_2 = 1 $, the dominant peaks and their ratios are found as $ 48/34 \approx 1.41 $ ($ \approx 3:2 $) and $ 99/48 \approx 2.06 $ ($ \approx 2:1 $). At $ S_1 = 2 $, $ S_2 = 1 $, the numerically obtained ratios are $ 32/23 \approx 1.39 $ ($ \approx 4:3 $) and $ 47/23 \approx 2.04 $ ($ \approx 2:1 $), indicating both resonant and overtone-type harmonics. In the case of $ S_1 = 4 $, $ S_2 = 1 $, the ratios $ 38/21 \approx 1.81 $ ($ 9:5 \approx 2:1 $) and $ 55.5/38 \approx 1.46 $ ($ \approx 3:2 $) are obtained. Finally, for $ S_1 = 5 $, $ S_2 = 1 $, the observed ratios are $ 33/26 \approx 1.27 $ ($ \approx 4:3 $) and $ 72/51 \approx 1.41 $ ($ \approx 3:2 $). These numerically observed mode pairs are in strong agreement with the 3:2 HFQPO pairs detected in BH X-ray binaries such as $ GRO J1655-40 $ ($ 300 / 450 $), $ XTE J1550-564 $ ($ 184/ 276 $), and $ GRS 1915+105 $ ($ 41/ 67 $) [68−71]. In addition, the $ 2:1 $ and $ 4:3 $ branches found in numerical simulations correspond well to transient resonant behaviors observed in systems such as $ XTE J1859+226 $ [72] and $ H1743-322 $ [73].

Compared to the Schwarzschild case, the Lee-Wick models produce a broader and richer PSD structure. Around the Schwarzschild BH, the QPO frequencies exhibit relatively low amplitudes dominated by LFQPOs. In contrast, the Lee-Wick parameters $ S_{1} $ and $ S_{2} $ reshape the spacetime geometry, amplifying the turbulence within the shock cone and enhancing the observability of strong, high-amplitude QPO peaks. Increasing $ S_{1} $ strengthens the gravitational confinement and slightly lowers the oscillation frequencies, whereas increasing $ S_{2} $ intensifies the oscillatory curvature component, generating more complex modulations and noticeable frequency upshifts. These results demonstrate that the Lee-Wick corrections produce observable differences in the PSD signatures, offering a promising avenue for identifying potential quantum-gravity imprints in accreting BH systems.

In the Block-2 configurations presented in Table 1, the Lee-Wick BH exhibits strong spacetime corrections. In this regime, both the exponential term governed by $ S_{1} $ and the oscillatory term determined by $ S_{2} $ become significantly more pronounced in shaping the metric function. The parameters $ S_{1} $ and $ S_{2} $ attain values that amplify the higher contributions from the derivative, thereby leading to substantial deviations from the Schwarzschild solution and, in some cases, to the formation of multiple horizons. The oscillatory component associated with $ S_{2} $ introduces spatially varying corrections to the gravitational potential, reflecting the characteristic behavior of complex mass poles in the Lee-Wick formulation. As a result, the curvature profile of spacetime may develop nontrivial features that encode the influence of quantum-scale modifications of gravity. The parameter values used in Block-2, therefore, enable us to explore how the fundamental structure of spacetime is altered in the strong Lee-Wick regime. As discussed in detail in the following, the numerical modeling of this regime reveals that strong Lee-Wick corrections modify the overall morphology of the accretion flow, alter the formation of the bow shock, and significantly influence the QPOs that arise within or near the shock region.

In order to reveal the effect of the parameter $ S_{2} $, which represents the oscillatory component of $ f(r) $, on the shock cone structure formed around the non-rotating Lee-Wick BH, the numerical results presented in Fig. 15 have been analyzed. In Fig. 15, using the initial conditions given in Block-2 of Table 1, the rest mass density color and contour plots, together with the velocity vector field, are presented to illustrate the morphology of the shock cone produced by BHL accretion and to demonstrate the influence of parameters $ S_{1} $ and $ S_{2} $ on the cone structure. In the upper-left snapshot, corresponding to $ S_{1}=1.5 $ and $ S_{2}=1.5 $, it is observed that both the exponential ($ S_{1} $) and oscillatory ($ S_{2} $) terms moderately affect the morphology of the shock structure. In this case, as a result of BHL accretion on the downstream side, a well-defined and symmetric shock cone similar to those of Block-1 (Fig. 9) is formed. Since $ S_{1} $ is relatively large, the exponential term becomes dominant and thus the Lee-Wick metric behaves similarly to the Schwarzschild case, producing a stable and narrow shock cone. In this regime, the spacetime curvature acts as a strong focusing potential with little oscillatory distortion.

Figure 15.  (color online) Same as Fig. 9, but this time the rest-mass density is shown, particularly near the BH, for the models in Block-2 listed in Tab. I. It can be seen that, depending on the model, the shock cone forms around the BH and consequently produces a bow shock or strong instabilities.

In the upper-right panel of Fig. 15, the cases $ S_{1}=1.5 $ and $ S_{2}=4 $ are modeled. Due to the larger value of $ S_{2} $, this configuration exhibits rapid fluctuations in curvature. As seen in Fig. 15, the cone becomes wider and encloses the BH, forming a bow-shock-like structure. The resulting pattern shows asymmetric density contours and visible instabilities near the cone’s apex. This behavior arises because the oscillatory curvature term governed by $ S_{2} $ generates varying local gravitational gradients, leading the flow to overshoot and produce quasi-periodic pressure waves.

In the lower-left panel of Fig. 15, the cases $ S_{1}=1 $ and $ S_{2}=1.5 $ are considered, where the lower value of $ S_{1} $ corresponds to a weaker exponential damping. Consequently, the Lee-Wick corrections remain strong even far from the horizon. The cone opening angle increases, and the density gradient becomes flatter, indicating a weakening of the gravitational confinement. However, since $ S_{2} $ is moderate, the oscillations remain mild and the flow is mostly stable, representing a transitional regime between Schwarzschild-like and oscillatory dynamics. In this case, a region of enhanced turbulence appears near the BH horizon.

In the lower-right panel of Fig. 15, the numerical results for $ S_{1}=1 $ and $ S_{2}=2 $, a configuration that favors strong Lee-Wick effects, are shown. Here, the exponential suppression is weak while the oscillatory curvature is prominent. This setup produces a dense, thick shock region that almost completely surrounds the BH, forming a bow shock. The formation of this bow shock arises from the collision of supersonically flowing matter with the high-pressure cavity in the downstream region. For larger $ S_{2} $ values, the enhanced oscillatory behavior modifies the effective potential around the Lee-Wick BH, increasing local compressions. This leads to nonlinear coupling between radial and azimuthal modes within the bow shock, giving rise to strong quasi-periodic oscillatory behavior. The development of these shock cones and bow shocks is in complete agreement with the theoretical analyses discussed above and with the effective potential, force profile, and ISCOs results presented in Fig. 2, Fig. 3, Fig. 4, and Fig. 5.

Analyzing the azimuthal variation of the rest-mass density is crucial for understanding the dynamical structure of the physical system forming in the strong gravitational field and near the ISCOs region, as well as for revealing the influence of the Lee-Wick parameters on spacetime geometry. Therefore, in Fig. 16, the azimuthal variation of the rest-mass density is presented. The top panel shows the behavior in the strong gravitational field at $ r = 3.56M $, while the bottom panel illustrates the distribution near the ISCOs at $ r = 6.11M $. In each panel, using the Block-2 initial conditions listed in Table 1, the wild behaviors arising from the spacetime modifications induced by the Lee-Wick parameters are demonstrated. As observed at $ r = 3.56M $, almost all models exhibit significant deviations compared to the Schwarzschild case. In particular, in the two models with higher values of oscillatory parameters ($ S_{2} = 2 $ and $ S_{2} = 4 $), pronounced and asymmetric peaks in density appear, revealing the formation of strong shock waves due to enhanced curvature oscillations. These strong fluctuations cause the accreting matter near the BH horizon to compress sharply, forming a bow-shaped shock cone surrounding the BH. In contrast, in cases with moderate or smaller $ S_{2} $ values ($ S_{2} = 1.5 $), the density profile becomes smoother, and a weaker shock wave is produced. This behavior originates from the exponential damping determined by $ S_{1} $, which leads to a more stable and symmetric accumulation of matter around the BH.

Figure 16.  (color online) Similar to Fig. 10, the azimuthal variation of the rest-mass density at two different radial locations is shown for the Block-2 models listed in Tab. I. In the upper panel, the variation of the rest-mass density at $ r = 3.56M $ is presented for different $ S_1 $ and $ S_2 $ values and compared with the Schwarzschild case. In the lower panel, the azimuthal behavior of the rest-mass density around the ISCOs, at $ r = 6.11M $, is displayed. It is clearly observed that, especially at $ r = 3.56M $, the region closest to the horizon, the behavior of the rest-mass density in Lee-Wick gravity significantly deviates from that of the Schwarzschild spacetime.

As seen in the lower panel of Fig. 16 at $ r = 6.11M $, corresponding to the ISCOs region, the influence of the Lee-Wick corrections is significantly weakened. Consequently, almost all models converge toward the Schwarzschild solution, although the maximum value of the rest-mass density remains relatively lower compared to the Schwarzschild case. However, the model with the highest oscillatory parameter value ($ S_{1}=1 $, $ S_{2}=2 $) still shows a clear deviation from Schwarzschild behavior. This indicates that increasing $ S_{2} $ not only alters the near-horizon structure but also significantly modifies the morphology of the physical flow even in slightly distant regions. These numerical results are in complete agreement with the theoretical predictions discussed earlier in the paper. Increasing $ S_{2} $ enhances the oscillatory mode of curvature and modifies the effective potential, resulting in stronger local gravitational gradients, nonlinear coupling between modes, and the formation of intense shock structures near the BH.

For the Block-2 models given in Table 1, the time evolution of the mass accretion rate calculated at $ r = 3M $ is presented in Fig. 17. As clearly seen in each model, different levels of instabilities are produced. The most stable configuration, characterized by the smallest oscillation amplitude, corresponds to the case with the lowest $ S_{2} $ value and the strongest exponential damping term. Its fluctuations are low in amplitude and quasi-regular, consistent with a spacetime metric closer to the Schwarzschild geometry and with weaker post-shock turbulence. In Fig. 17, when $ S_{1} $ is fixed and $ S_{2} $ increases, the signal becomes noticeably noisier and broadband, with significantly larger oscillation amplitudes compared to the other models. The enhanced oscillatory component of the lapse function $ f(r) $ drives stronger curvature ripples, which feed variability into the shock cone/bow-shock system and the downstream cavity. This leads to the excitation of strong global oscillation modes.

Figure 17.  (color online) For the models given in Block-2 of Table 1, the time evolution of the mass accretion rate is presented. The figure illustrates how the oscillation strength of the shock cone and bow shock varies with different values of the Lee-Wick gravity parameters $S_{1}$ and $S_{2}$. Increasing the parameter $S_{2}$ enhances the oscillatory behavior of the accretion flow, leading to stronger temporal fluctuations and quasi-periodic variations in the mass accretion rate.

When the exponential term is suppressed (i.e., smaller $ S_{1} $), both the mean accretion rate and its temporal variation increase since the effective gravitational force becomes more attractive and the infalling matter is compressed more efficiently before forming a shock. As a result, larger values of $ S_{2} $ amplify the oscillatory nature of the metric and destabilize the post-shock flow. At the same time, as shown in the theoretical analyses presented earlier, increasing $ S_{1} $ or $ S_{2} $ deepens the effective potential, thereby enhancing the gravitational attraction. As also predicted by theory, stronger compression near the BH leads to shifts in the ISCOs position, altogether yielding the stronger, less-regular variability observed in the high $ S_{2} $ runs.

The oscillatory curvature creates rapidly varying local gravitational gradients and, as a result, strong instabilities appear, as evident in Fig. 17. The accretion flow repeatedly overshoots, reflects within the high-pressure cavity, and drives standing and advected shock waves, leading to large-amplitude and irregular fluctuations in the mass accretion rate. This mechanism also explains the different shock morphologies observed in Fig. 15. In all model configurations, the generated instabilities represent key features that support the formation of QPOs. In particular, the strong oscillations observed in the high $ S_{2} $ cases correspond to quasi-periodic variations, confirming the emergence of nonlinear oscillatory behavior consistent with the QPO formation scenario discussed in Section V B 2.

As seen in Fig. 17, the mass accretion rate exhibits strong oscillatory behavior depending on the parameters $ S_{1} $ and $ S_{2} $, reflecting the dynamical influence of the shock cone or the bow-shock structure formed around the BH. For each Block-2 model with different $ S_{1} $ and $ S_{2} $ values, PSD analyses are performed using the mass accretion rate data shown in Fig. 17. As clearly observed, each configuration displays a rich harmonic structure, and the resulting QPO analyses are presented in Fig. 18. In Fig. 18, for the cases of $ S_{1}=1.5 $ and $ S_{2}=1.5 $, the PSD analysis reveals that the QPO frequencies appear in the range of $ 7-70 $ Hz. The frequency ratios are approximately $ 2:1 $ $ \approx 19/10 $, $ 3:1 $ $ \approx 30.6/10 $, and $ 4:1 $ $ \approx 36/10 $, corresponding to weakly coupled LFQPOs in a spacetime similar to Schwarzschild. When $ S_{2} $ increases to $ 4 $, the dominant frequency shifts to about $ 12 $ Hz, with several overtones that extend to $ \sim70 $ Hz. The characteristic ratios in this model are $ 3:2 $ $ \approx 17.8/12.4 $, $ 3:1 $ $ \approx 36.5/12.4 $, $ 4:1 $ $ \approx 51.4/12.4 $, and $ 6:1 $ $ \approx 73.8/12.4 $, indicating the presence of a strong nonlinear coupling between oscillation modes.

Figure 18.  (color online) PSD analysis of the mass accretion rate for the initial conditions given in the Block-2 configuration shows the dependence on the Lee-Wick parameters. For each model, a rich set of oscillation frequencies is observed, clearly demonstrating the influence of the Lee-Wick parameters on the resulting QPO frequencies. Increasing the parameter $S_{2}$ shifts the dominant frequencies toward the HFQPO regime, indicating stronger nonlinear coupling between the radial and azimuthal oscillation modes within the accretion disk.

For the configurations $ S_{1}=1 $ and $ S_{2}=1.5 $, the PSD analysis shows slightly higher fundamental modes in the range of $ 11-20 $ Hz, suggesting a moderate oscillatory modulation of the accretion flow. The observed frequency ratios are $ 3:2 $ $ \approx 20.7/13.6 $, $ 2:1 $ $ \approx 25.8/13.6 $, and $ 5:2 $ $ \approx 33.8/13.6 $, revealing harmonically related overtones. However, the most striking behavior appears in the cases of $ S_{1}=1 $ and $ S_{2}=2 $, where the PSD displays peaks extending up to nearly $ 800 $ Hz, with dominant features at $ 86 $, $ 136 $, $ 223 $, $ 285 $, and $ 324 $ Hz. These yield characteristic ratios of $ 3:2 $, $ 5:2 $, $ 3:1 $, and $ 4:1 $, which are typical of HFQPO.

These results clearly show that the parameter $ S_{1} $ controls the exponential damping of the quantum corrections, enhancing the stability of accretion and favoring the appearance of low-frequency modes. In contrast, increasing $ S_{2} $ strengthens the oscillatory component of the metric, amplifies nonlinear curvature coupling, and shifts the resulting oscillation frequencies toward the HFQPO regime. The transition from LFQPO to HFQPO in the simulations corresponds to the inward shift of the effective potential and ISCOs as $ S_{2} $ increases, leading to faster epicyclic oscillations consistent with theoretical predictions of the Lee-Wick spacetime.

Finally, the QPOs obtained from the numerical simulations confirm that the Lee-Wick parameters naturally allow for the coexistence of multiple QPO regimes. The LFQPOs observed in low-$ S_{2} $ cases, typically between $ 7 $ and $ 30 $ Hz, are comparable to those detected in black-hole X-ray binary systems. In contrast, high-$ S_{2} $ models produce QPOs in the range of $ 80 $ to $ 300 $ Hz, consistent with observed HFQPOs, such as $ \sim67 $ Hz in GRS 1915+105, $ \sim184/276 $ Hz and $ \sim300/450 $ Hz $ 3:2 $ pairs in GRO J1655-40 and XTE J1550-564 [68−71]. In particular, the $ 136/86 \approx 3:2 $ ratio obtained for the $ S_{1}=1 $, $ S_{2}=2 $ model closely matches the observational data from GRS 1915+105, confirming the physical relevance of the simulated QPO frequencies to real astrophysical systems.

According to the numerical results presented in Section V, the morphology of the accretion flow and its time-dependent variations around the nonrotating Lee-Wick BH are strongly dependent on the two coupling parameters, $S_{1}$ and $S_{2}$, which characterize the space-time metric of this BH. While $S_{1}$ represents the exponential component of the metric function $f(r)$, $S_{2}$ corresponds to its oscillatory component. In the parameter region defined as Block-1 in Table 1, where the corrections are weak ($S_{1} \geq 2$ and $S_{2} \leq 1.5$), the exponential damping dominates, causing the spacetime geometry to behave almost identically to that of the Schwarzschild solution. Consequently, the BHL accretion pattern forms a smooth, narrow, and quasi-steady-state shock structure. The corresponding temporal evolution of the mass-accretion rate exhibits mild periodic modulations with low amplitude, producing LFQPOs-like features with harmonics around a few tens of hertz. These results are consistent with type-C LFQPOs observed in Galactic black-hole binaries such as XTE J1550-564 and GRS 1915+105. Therefore, it can be inferred that the QPO signals detected from these sources likely originate from stable shock oscillations and weak nonlinear coupling effects within the accretion flow in a spacetime where the Lee-Wick corrections are negligible.

In contrast, when the Lee-Wick metric is characterized by lower values of $S_{1}$ and higher values of $S_{2}$, strong oscillations arise, leading to rapid curvature fluctuations and even reversals of the gravitational gradient. As a result, the accretion flow becomes highly dynamical, forming broad and asymmetric shock cones as well as distinct bow-shock morphologies. The generated bow shock entirely encloses the BH, enhancing the strength of the oscillatory behavior. The temporal evolution of the accretion rate in this regime reveals strong, short-time-scale quasi-periodic modulations. The PSD analysis demonstrates that dominant peaks occur within the $ 10-800 $ Hz range. These peaks exhibit frequency ratios such as $ 3:2 $, $ 2:1 $, $ 3:1 $, and $ 5:2 $, which arise due to the strong nonlinear coupling between the radial and azimuthal oscillation modes in the region of intense dynamical activity, namely, within the shock cone and bow-shock zones. The numerical results clearly indicate that as $S_{2}$ increases, the characteristic QPOs systematically shift toward the HFQPO regime. These findings are in full agreement with the theoretical results presented in Section IV, where it was shown that the oscillatory curvature enhances the radial epicyclic frequency and causes the ISCOs to move inward. In physical terms, a larger $S_{2}$ increases the oscillatory component of the gravitational potential, thereby amplifying the nonlinear mode coupling and pushing the disk oscillations into the relativistic regime.

When the agreement between the numerical results and observations is considered, the frequency ranges obtained from the simulations show strong consistency with those observed in microquasars and X-ray binary systems. In the Block-1 regime, the LFQPOs identified in the $ 5-30 $ Hz range are consistent with the type-B/C QPOs observed in sources such as XTE J1550-564 [78], GX 339-4 [77], and GRS 1915+105 [76]. In contrast, HFQPOs found in the Block-2 regime, within the $ 80-300 $ Hz range, coincide with the $ 67 $ Hz oscillation of GRS 1915+105, the $ 184/276 $ Hz pair of XTE J1550-564 [75], and the pair of $ 300/450 $ Hz detected in GRO J1655-40 [74]. In particular, the $ 3:2 $ frequency ratio observed in the PSD analysis for the cases $S_{1}=1$ and $S_{2}=2$ exhibits a strong parallel to the resonance ratios reported in these observational sources. Moreover, the amplitude growth and spectral broadening observed in the numerical simulations for large $S_{2}$ values are remarkably similar to the irregular variability and turbulence-induced transitions detected during state changes of transient BH binary systems.

The shift of the numerically extracted QPOs can also be understood from the theoretical analysis of the effective potential and ISCOs structure presented in Section III. As shown in the paper, increasing $S_{1}$ deepens the minimum of the potential through exponential suppression, while increasing $S_{2}$ superposes an oscillatory term on the potential. This addition modulates the curvature and splits the potential well into sub-minima. The strong interplay between these two terms naturally leads to the formation of multiple resonance cavities, within which matter undergoes alternating phases of compression and rarefaction. As a result, both LFQPOs and HFQPOs are effectively excited. This behavior closely parallels discrete peaks and broadband noise observed in systems such as H1743-322 [79] and 4U 1630-47 [80], where the presence of multiple QPO peaks has been confirmed observationally. Therefore, the numerical simulations performed in this study provide a physically motivated mechanism through which higher-derivative gravitational corrections, parameterized by $(S_{1}, S_{2})$, can leave measurable imprints on the temporal spectra of accreting BHs.

From the results obtained in the numerical simulations, two distinct observational implications can be identified. First, according to the results derived within the Lee-Wick framework, the QPO frequencies formed around nonrotating BHs depend not only on the mass of the BH but also on the parameters $S_{1}$ and $S_{2}$. This introduces new degrees of freedom that could be constrained by future X-ray timing observations. Next-generation missions such as eXTP [81], STROBE-X [82], and Athena [83] may be able to test whether the observed frequency ratios and their energy dependencies correspond to the parameter dependence on $S_{1}$ and $S_{2}$ discussed in detail here. Second, the morphological structure of the bow shock that emerges at large $S_{2}$ values may manifest itself as asymmetric brightness features or variable absorption patterns in Event Horizon Telescope scale [84] images of accretion flows, offering a potential multiwavelength diagnostic of quantum-gravity corrections. Furthermore, the oscillatory curvature and enhanced compression that develop particularly near the horizon could affect the jet-launching conditions, potentially linking the Lee-Wick parameter $S_{2}$ to episodic radio flaring similar to that observed in GRS 1915+105 [85]. Lastly, the Lee-Wick gravity framework should be viewed not merely as a theoretical model capable of explaining QPO formation but also as a promising observational channel through which the effects of higher-derivative quantum gravity could be tested in future astrophysical observations.

In the weak Lee-Wick region, defined by $ S_{1} \ge 2 $ and $ S_{2} \le 1.5 $, the exponential suppression term in the metric function $ f(r) $ dominates over the oscillatory correction term. The reasoning behind the sign choices for the parameters $ S_1 $ and $ S_2 $ has been discussed [21]. In this work, we have deliberately explored parameter values both inside the admissible region defined by the constraint $ S_2 \gt 0 $ and $ -2\sqrt{S_2} \lt S_1 \lt 2\sqrt{S_2} $ and beyond the constraint, to sharply reveal how these two parameters shape our analysis. Interestingly, we have found both kinds of values physically viable as beyond the constraint we obtained almost similar results to those observed with the constraint $ S_2 \gt 0 $ and $ -2\sqrt{S_2} \lt S_1 \lt 2\sqrt{S_2} $. There was only a bigger difference in the case of shock cone analysis. As a consequence, the numerical simulations demonstrate that the shock cone morphology produced in these models does not differ significantly from that of the Schwarzschild spacetime. The curvature of spacetime remains regular and exhibits no substantial deviation from general relativity, leading to a smooth and well-defined accretion structure. As in the standard Schwarzschild BH, the BHL mechanism produces a symmetric shock cone in the downstream region, which maintains a quasi-steady configuration over time while undergoing mild oscillations. These oscillations originate from small perturbations within the post-shock region and appear as coherent harmonic variations in the mass-inflow rate. Furthermore, the opening angle of the shock cone is found to be in close agreement with the results obtained for the Schwarzschild case. This confirms that the Lee-Wick corrections in this regime are weak, the influence of the higher-derivative terms is minimal, and the overall accretion dynamics and the formation of the shock cone are entirely governed by standard general relativistic gravity.

A detailed analysis of test particle motion around the non-rotating Lee-Wick BH shows that parameters $ S_1 $ and $ S_2 $ are fundamental in determining the spacetime orbital structure. The effective potential shows clear distinctions between stable and unstable circular orbits and provides a consistent criterion for orbital stability through its minima and maxima. The profiles of energy and angular momentum show that for nearly all parameter regimes, $ S_1 $ and $ S_2 $ can be increased, and the energy needed for circular motion is lowered. The same modest impact (decrease) is seen for angular momentum, particularly in certain parameter regimes. These tendencies are seen in the radii of the ISCOs, which for large $ S_1 $ and $ S_2 $ values show the same trend of decreasing close to a saturation value. In all cases, the Lee-Wick modifications exhibit clear and meaningful fundamental changes to the stability and energy criterion and radius of characteristic orbits. Overall, the underlying quantum-inspired corrections to geodesic motion in this BH spacetime are underscored by the clear, fundamental, and physical changes to the geodesic motion.

In the strong Lee-Wick region, where $ S_{1} $ is small and $ S_{2} $ is large, the exponential and oscillatory components of the metric function contribute with comparable strength, leading to pronounced modifications in the curvature structure of spacetime. As revealed by the numerical simulations, these strong gravitational corrections produce rapid oscillations in curvature and abrupt variations in the metric gradient near the BH horizon. Such variations strongly influence the dynamics of the accreting matter, causing the downstream shock cone to widen, become asymmetric, and exhibit an increasing level of instability over time. Eventually, the initially well-defined shock cone evolves into a bow shock-like configuration that fully surrounds the BH, signifying the transition to a qualitatively different physical regime dominated by oscillatory spacetime curvature. This behavior clearly demonstrates that the results in the strong Lee-Wick domain deviate substantially from the Schwarzschild case, as the accretion flow becomes turbulent and the governing dynamics are dominated by nonlinear gravitational couplings. In particular, the time-dependent evolution of the mass accretion rate in this region shows highly irregular oscillations occurring on very short timescales. The PSD analyses reveal that these intense instabilities generate multiple frequency components in the range of $ 10 $-$ 800\,\mathrm{Hz} $, producing a spectrum of QPOs. These signatures confirm that strong mode coupling arises due to the oscillatory curvature term governed by $ S_{2} $, resulting in QPOs whose amplitudes and frequency distributions depend sensitively on the Lee-Wick parameters. Therefore, this strong-correction regime provides an important theoretical framework for probing the influence of quantum-gravity effects on the spacetime geometry around BHs, offering a promising avenue for testing the astrophysical implications of Lee-Wick gravity.

The numerical results have shown that the astrophysical manifestations of Lee-Wick gravity vary strongly depending on the parameters $ S_{1} $ and $ S_{2} $. In the weak correction limit, the morphology of the shock cone closely resembles that of the Schwarzschild solution. In this case, the mass accretion flow exhibits smooth behavior with mild temporal modulations capable of explaining the LFQPOs. The numerical findings are in good agreement with the QPOs observed in X-ray binary systems such as XTE J1550-564 and GRS 1915+105. The presence of such LFQPOs indicates that small perturbations within the post-shock region are coherently sustained as harmonic oscillations of the accretion flow, whereas the influence of higher-derivative corrections remains minimal. Therefore, the weak field limit can account for observed astrophysical systems that exhibit stable, low-amplitude QPOs. Furthermore, the stability and regularity of the shock cone in this regime suggest that the observed frequency distribution and amplitude of LFQPOs may impose constraints on the Lee-Wick parameters, potentially setting limits on $ S_{1} $ and $ S_{2} $ for such sources. However, the numerical results obtained for the strong correction regime reveal distinct observational consequences. The interaction between the spacetime parameters in this domain, characterized by nonlinear gravitational couplings and strong curvature oscillations near the BH horizon, can explain the HFQPOs observed in various sources. In this region, the intense and rapidly varying oscillations provide the physical conditions necessary for HFQPO formation. The PSD analyses of the accretion rate demonstrate that the QPO frequencies extend from a few tens to several hundred hertz. Numerically, these HFQPOs arise from the coupling between radial and azimuthal modes, producing frequency ratios such as $ 3:2 $, $ 2:1 $, and $ 5:2 $. These ratios are in excellent agreement with those observed in GRO J1655-40 (300/450\,Hz) and XTE J1550-564 (184/276\,Hz). Hence, the oscillatory curvature term governed by $ S_{2} $ drives the emergence of nonlinear resonant modes, leading to the transition from LFQPOs to HFQPOs. Consequently, while the weak Lee-Wick regime reproduces the classical LFQPO phenomenology, the strong correction domain naturally explains the appearance of HFQPOs and broadband oscillatory structures. Thus, within a single theoretical framework, Lee-Wick gravity can unify both LFQPO and HFQPO behaviors, offering a valuable opportunity to investigate quantum gravity-induced spacetime modifications in the vicinity of BHs.