An Advanced Fuel Efficiency Optimization Model with Fractional Programming

Md Isfakul Anam,  Clarkson University, Tuyen Vu,  Clarkson University

Md Isfakul Anam, T. T. Nguyen, and T. Vu with Clarkson University, Potsdam, NY, USA; Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding Author: T. Vu, Email<EMAIL_ADDRESS>
Reducing the fuel consumption within a power network is crucial to enhance the overall system efficiency and minimize operating costs. Fuel consumption minimization can be achieved through different optimization techniques where the output power of the generators is regulated based on their individual efficiency characteristics. Existing studies primarily focus either on maximizing the efficiency function or minimizing the operating cost function of the generators to minimize fuel consumption. However, for practical implementation, it becomes imperative to incorporate a function within the optimization framework to represent the fuel consumption rate directly. This study introduces a novel approach by formulating a minimization problem with a sum-of-ratios objective function representing the fuel consumption rate. However, optimization problems with sum-of-ratios objective functions or constraints are extremely challenging to solve because of their strong nonlinearity. To efficiently solve the formulated problem, a fractional programming (FP) approach is adopted in this study. This reformulation technique significantly reduces the solution time of the optimization problem and provides a better solution than nonlinear programming (NLP). In addition, the reformulated problem can also be applied to large-scale systems where the NLP fails to converge. The proposed methodology of this study is tested on the notional MVAC ship system, modified IEEE 30-bus and IEEE 118-bus systems. The results demonstrate that the model successfully minimizes fuel consumption by effectively scheduling the generator and ESS dispatch.

Fuel consumption, system efficiency, energy management system, sum-of-ratio problem, fractional programming.

[01]\(P_{i,g}, Q_{i,g}\)Real and reactive power output of $i$-th generator
[02]\(P_{i,l}, Q_{i,l}\)Real and reactive load at $i$-th bus
[03]\(P_{i,inj}, Q_{i,inj}\)Real and reactive power injection at $i$-th bus
[04]\(V_i, \theta_{ik}\)voltage and voltage angle difference
[05]\(G_{ik}, B_{ik}\)Conductance and susceptance of line $ik$
[06]\(P_{ik}\)Real power capacity of line $ik$
[07]\(DR_i, UR_i\)Down and up ramp rate of $i$-th generator
[08]\(P_i^{C,max}, P_i^{D,max}\)Maximum charging and discharging rate of $i$-th ESS
[09]\(E_{i,b}\)Energy stored at $i$-th ESS
[10]\(P_{i,b}^r\)Output power of $i$-th ESS
[11]\(SOC_i\)State of charge of $i$-th ESS
[12]\(\alpha\)fuel energy density
[13]\(N_g\)Number of generator
[14]\(T\)Total planning horizon
[15]\(\triangle t\)Time step
[16]\(B\)Number of buses

§ INTRODUCTION

Continuous increments in load demand on power systems due to the growing number of consumers impose a great challenge for the engineers and operators in the field. It is essential to supply the increased load requirements of a power system territory by introducing additional energy sources and/or expanding the capacity of the existing sources [1]. Both solutions result in higher fuel costs incurred for operating a large number of distributed generators. Furthermore, the survival period of stand-alone power networks, such as islanded microgrids or ship electrical power systems, largely depends on the fuel consumption of the distributed generators. As a result, a state-of-the-art system efficiency model is required to optimize the fuel consumption of a system.

System efficiency optimization techniques refer to minimizing the fuel consumption rate by regulating the output power of the generators within the system. Over the past decades, a significant number of research has been conducted to improve system efficiency using various optimization methods.
One notable contribution is presented in [2], where the authors propose a novel approach to analyze the characteristics of the efficiency function, leading to the determination of the maximum total power supply and overall efficiency for such systems. In [3], a genetic algorithm modeling framework is presented where the optimization problem involves minimizing the thermal cost function. The cost function is constructed based on the power generation characteristics of the hydropower plant. A similar approach can be found in [4], where the daily optimal generation scheduling problem (DOHGSB) is solved by implementing a unique differential evolution algorithm. The objective function is formulated by analyzing the hydropower plant characteristics or input-output curve.
In [5], the authors introduce a novel distributed algorithm to maximize system efficiency. A fourth-order and a third-order efficiency function for the main and auxiliary power generation module (PGM), respectively, are optimized using a distributed crow search algorithm (DCSA).
The approaches in [2]-[5] focus on optimizing the efficiency function or the generation characteristics of the generators, which can be complex to apply to systems consisting of generators with different ratings. Also, since the generator efficiency function is not a straightforward representation of the system's fuel consumption, it is impractical to optimize the efficiency function with an objective to improve fuel efficiency.

Some literature can be found where the optimization problem is formulated to minimize the fuel consumption rate or fuel cost directly. In [6], four power-sharing schemes are presented to establish a unit commitment strategy to minimize fuel costs. However, instead of considering the generator efficiency function, the authors utilized a typical fuel consumption characteristic curve, which largely depends on the capacity and model of a generator.
The authors in [7] implement a recursive method to estimate a second-order polynomial model of specific fuel consumption. The model is later used to determine the optimal load distribution between the different generators. A similar approach is found in [8], where dynamic programming is used to solve the formulated problem. However, both optimization models are suitable for simple power networks since they don't include AC power flow or energy storage models.
A minimum hourly fuel consumption curve interpolated by a quadratic equation is used in [9] to minimize the fuel consumption of hybrid electric vehicles. This method can not be extended for transmission or distribution systems since most power system constraints are not included in the model.
In [10], the authors utilize a 3D map of brake-specific fuel consumption (BSFC) in terms of the rotating speed of the drive and generated mechanical torque to determine the minimum point of diesel engine (DE) fuel consumption. Authors in [11] also apply a similar method where a speed vs. power curve of the diesel engine (DE) is exploited to achieve minimum fuel consumption conditions. Nonetheless, these methods have two major drawbacks: they are implemented particularly for the DC ship systems, and the strong nonlinearity of the BSFC curve increases the complexity of determining the optimal operation.

Reduction in fuel cost can also be achieved indirectly through economic dispatch (ED) optimization problems [12], [13]. In the ED problem, a polynomial objective function (generally in quadratic form) representing the cost of the generator dispatch is optimized to supply the demand most economically. However, since the objective function of the ED problem does not represent the fuel consumption rate, these formulations are unable to accomplish the highest system efficiency. They are also inconvenient for long-term planning of fuel usage and impractical to apply to systems where achieving optimal fuel consumption rate is the main goal.

Although the strategies discussed above for improving overall system performance and reducing fuel cost have their own merits, the key to minimizing fuel consumption lies in introducing a function that accurately represents the fuel consumption rate. By directly representing fuel consumption in the optimization process, researchers can effectively tackle the core challenge of reducing fuel usage and achieving greater energy efficiency in the system. Therefore, future studies may benefit from exploring methodologies considering this crucial aspect when addressing system efficiency optimization.

This paper presents a novel approach by introducing a unique sum-of-ratios objective function, which directly represents the fuel consumption rate. Unlike conventional methods that optimize the polynomial generator efficiency function or cost of operation, this sum-of-ratios formulation offers a practical and more efficient solution. However, solving multiple ratio optimization problems, known as Fractional Programming (FP), has been proven to be NP-hard [14]. The convergence of these problems with established nonlinear optimization methods can take an extensive amount of time. In addition, when the sum-of-ratios problem involves more than 20 ratios, the current approaches struggle to find a solution within a reasonable timeframe [15],[16].
Due to these challenges, directly solving the sum-of-ratios optimization problem for energy management systems (EMS) is not feasible. Especially this method will be inapplicable for large-scale systems with a high number of variables. It necessitates the development of an efficient reformulation and solution technique to address multiple ratio problems effectively. Finding innovative approaches to tackle these difficulties will be crucial in making the proposed sum-of-ratios method a practical and scalable solution for optimizing fuel consumption in real-world energy systems.

A substantial body of literature exists on Fractional Programming (FP), but the emphasis has primarily been on single-ratio problems. Among the renowned reformulation techniques, the Charnes-Cooper Transform, [17] [18] is notable for proposing an algorithm to solve single ratio linear FP problems by introducing two new variables and converting the fractional problem into a linear problem. Another classical technique, Dinkelbach's Transform [19], reformulates the single ratio problem using a new auxiliary variable updated iteratively until convergence is achieved.
[20] presents a formulation of a linear problem equivalent to a single ratio linear FP problem where some duality properties are used to prove the equivalence. For quadratic FP problems, where both the numerator and denominator are quadratic functions, a new method called the decomposition fractional separable method is proposed in [21] using linear programming techniques.
An alternative approach to solving single-ratio quadratic FP is outlined in [22], employing Taylor series expansion for effective reformulation.

The literature discussed in references [17]-[22] primarily focused on solving single-ratio FP problems and cannot be directly extended to handle multi-ratio problems. Addressing multiple ratio problems, as encountered in the sum-of-ratios function, remains a challenge that requires innovative and efficient reformulation and solution techniques. However, in [23], the authors proposed an extension of Dinkelbach's Transform specifically tailored to address multi-ratio FP problems. Nonetheless, this method was later refuted by Falk and Palocsay [24], who demonstrated its limitations through a numerical example.

To find the globally optimal solution for the sum-of-ratios problem, [25] introduced a practical method that involves solving a sequence of convex programming problems. In [26], a convexification strategy was employed to decompose fractional terms into convex and concave components. Then, a piecewise linearization technique was applied to approximate the concave terms effectively.
Additionally, [27] proposed a quadratic transform to tackle concave-convex multiple ratio minimization problems. In the case of generalized convex multiplicative functions, a reformulation technique was presented in [28], where the main problem was reformulated as a concave minimization problem with 2p variables. This reformulation technique could also be applied to sum-of-ratios FP problems if the multiplicative terms were replaced with a convex over a concave function [29].

In our study, the objective function is defined as a sum-of-ratios minimization problem with non-negative-convex numerator and positive-concave denominator terms. Due to this specific form of the problem, an appropriate algorithm should be selected to solve the formulated multiple ratios FP problem effectively. As a result, the reformulation technique presented in [28] is adopted in this paper. By leveraging this reformulation technique, the complexities of the sum-of-ratios problem can be effectively addressed, and an optimized solution can be found with a feasible convergence time.

The contributions of this paper are the following:

* A novel fractional objective function is introduced in this literature, which directly represents the fuel consumption rate of the generators. Unlike typical system efficiency optimization problems that use the efficiency function or the operating cost function as the objective, this unique formulation directly accounts for fuel consumption. This approach proves to be more efficient and practical compared to previous studies since it directly targets the core issue of minimizing fuel usage and improving overall system efficiency.

* To address the optimization problem in this study, the sum-of-ratios fractional programming (FP) algorithm is employed. to the best of our knowledge, this literature represents the first application of the FP method to solve the optimization problem for EMS efficiently. The reformulation technique with FP can also be applied to different power or communication system research where sum-of-ratios functions are used.

* The successful application of the FP algorithm, combined with the convex relaxation of nonlinear constraints, demonstrates that the proposed model is suitable for handling large-scale systems. This capability is exemplified through the model's effective implementation on the IEEE 118-bus system. By demonstrating its applicability to such a complex and extensive system, the paper establishes the scalability of the proposed approach for real-world energy management scenarios.

The remainder of the paper is organized as follows: the fuel efficiency problem formulation with sum-of-ratios objective function and its convex reformulation are presented in Section <ref>. In Section <ref>, the solution algorithm for the reformulated problem is described. The results for the notional MVAC ship system, IEEE 30-bus, and IEEE 118-bus system, accompanied with the performance comparisons, are demonstrated in Section <ref>. Finally, Section <ref> represents the conclusion and future work.

§ PROBLEM FORMULATION

§.§ Optimization Model for System Efficiency

This section presents a unique sum-of-ratios objective function for the optimization problem that directly represents the fuel consumption of the generators. The objective of the minimization problem is to minimize the fuel consumption rate of the generators over the planning horizon, which will maximize the system efficiency. The fuel consumed by a generator can be expressed by taking into account the generator's efficiency and the output power it produces. The objective function is the following:

\begin{equation}
\text{minimize } f = \frac{1}{\alpha}\sum_{t=0}^{T} \sum_{i\in N_g} \frac{P_{i,g}^t}{\eta_{i,g}} \Delta t.
\label{eq:System_Efficiency:Eq1}
\end{equation}

where, generator efficiency, $\eta = a_i p_{i,g}^2 + b_i p_{i,g} + c_i$ with $a$, $b$, and $c$ are generator specific constants, $\alpha$ is the fuel energy density (MWh/L),
$p_{i,g}$ is the per-unit output of the i-th generator: $p_{i,g} = {P_{i,g}}/{P_{i,b}}$;
$P_{i,g}^t$ is the generator output power at time t, and $P_{i,b}$ is the base power of i-th generator.
$N_{g}$ = total number of generators, $T$ = planning horizon, $\triangle t$ = each time period.
The following active and reactive power balance constraints are associated with the system:

\begin{equation}
P_{i,inj} = P_{i,g} - P_{i,l}
\label{eq:System_Efficiency:Eq2}
\end{equation}

\begin{equation}
Q_{i,inj} = Q_{i,g} - Q_{i,l}
\label{eq:System_Efficiency:Eq3}
\end{equation}

\begin{equation}
P_{i,inj} = \sum_{k\in B} V_i V_k (G_{ik} cos\theta_{ik} + B_{ik} sin\theta_{ik})
\label{eq:System_Efficiency:Eq4}
\end{equation}

\begin{equation}
Q_{i,inj} = \sum_{k\in B} V_i V_k (G_{ik} sin\theta_{ik} + B_{ik} cos\theta_{ik})
\label{eq:System_Efficiency:Eq5}
\end{equation}

where (<ref>) and (<ref>) are the AC power flow constraints for the system.
The following constraints should be included in the problem formulation to maintain the operational limits of the system:

\begin{equation}
P_{i,g}^{min} \leq P_{i,g}^t \leq P_{i,g}^{max},
\label{eq:System_Efficiency:Eq6}
\end{equation}

\begin{equation}
Q_{i,g}^{min} \leq Q_{i,g}^t \leq Q_{i,g}^{max},
\label{eq:System_Efficiency:Eq7}
\end{equation}

\begin{equation}
V_{i}^{min} \leq V_{i}^t \leq V_{i}^{max},
\label{eq:System_Efficiency:Eq8}
\end{equation}

\begin{equation}
\theta_{i}^{min} \leq \theta_{i}^t \leq \theta_{i}^{max},
\label{eq:System_Efficiency:Eq9}
\end{equation}

\begin{equation}
-P_{ik} \leq P_{ik}^t \leq P_{ik},
\label{eq:System_Efficiency:Eq10}
\end{equation}

\begin{equation}
-Q_{ik} \leq Q_{ik}^t \leq Q_{ik},
\label{eq:System_Efficiency:Eq10.1}
\end{equation}

\begin{equation}
-DR_{i} \leq P_{i,g}^{t+1} - P_{i,g}^{t} \leq UR_{i}.
\label{eq:System_Efficiency:Eq11}
\end{equation}

where (<ref>) and (<ref>) represent the generators' real and reactive power generation limits, (<ref>) and (<ref>) are the voltage and voltage angle limits, (<ref>) and (<ref>) are the line limits for real and reactive power, and (<ref>) is the ramp rate limit.

The energy storage system (ESS) plays a vital role in minimizing the fuel consumption by the generators. The following Energy Storage System (ESS) constraints are included in the optimization problem:

\begin{equation}
E_{i,b}^{t} = E_{i,b}^{t-1} - \eta_b P_{i,b}^{r,t} \triangle t
\label{eq:System_Efficiency:Eq12}
\end{equation}

\begin{equation}
-P_{i}^{C,max} \leq P_{i,b}^{r,t} \leq P_{i}^{D,max},
\label{eq:System_Efficiency:Eq13}
\end{equation}

\begin{equation}
\sum_{t=0}^{t} P_{i,b}^{r,t} = 0,
\label{eq:System_Efficiency:Eq14}
\end{equation}

\begin{equation}
SOC_{i}^{min} \leq SOC_{i}^{t} \leq SOC_{i}^{max},
\label{eq:System_Efficiency:Eq15}
\end{equation}

where (<ref>) indicates the energy conservation constraint, (<ref>) is the limit for charging or discharging rate, and (<ref>) is the state of charge (SOC) limit of the ESS. For the ESS, although the SOC can vary from 0 to 1 (0% to 100%), fully discharging can damage the battery permanently and shorten the life cycle of the battery [30]. In this paper, the minimum SOC is selected as 0.2 (20%). Eq.(<ref>) ensures that the sum of the total charging and discharging power over a planning period will be zero, which helps the system to recharge the battery before the next planning cycle.

§.§ Reformulation of the Objective Function

The objective function (<ref>) for the formulated optimization problem is a highly nonlinear sum-of-ratios case. In this section, the objective function is reformulated to a convex function using a technique used for solving generalized convex multiplicative problems. Later, the convex minimization method is utilized to solve the reformulated problem iteratively.

The general convex multiplicative minimization problem has the following structure:

\begin{equation}
\text{minimize } h(x) + \sum_{i=1}^{p} f_i(x)g_i(x)
\label{eq:System_Efficiency:Eq16}
\end{equation}

subject to  $x \in X$

where $h, f_i(x)$ and $g_i(x)$ for all $i$ are convex functions and $X\subset R^n$ is a convex set.
If $h(x)=0$, $f_i(x)=A_i(x)$, and $g_i(x)=1/B_i(x)$, (<ref>) will be in the following form:

\begin{equation}
\text{minimize } H(x) = \sum_{i=1}^{m} \frac{A_i(x)}{B_i(x)}
\label{eq:System_Efficiency:Eq17}
\end{equation}

subject to  $x \in X$
which is a sum-of-ratios problem, where $A_i(x)$ are non-negative, convex and $B_i(x)$ are positive, concave functions for all $i$.

The authors in [28] defined the following problem by introducing 2m auxiliary variables $\zeta_i$ and $\beta_i$, where $i=1,2,3,....m$:

\begin{equation}
\text{minimize } F(x, \zeta, \beta) = \frac{1}{2} \sum_{i=1}^{m} [\zeta_i (A_i(x))^2 + \beta_i (B_i(x))^2]
\label{eq:System_Efficiency:Eq18}
\end{equation}

subject to  $x \in X$
$\zeta_i \beta_i \ge 1$
$(\zeta, \beta) > 0$

where $\zeta = (\zeta_1, \zeta_2, ....\zeta_m)$ and $\beta = (\beta_1, \beta_2, ..... \beta_m)$. It can be proved that, if $(x^*, \zeta ^*, \beta ^*)$ is an optimal solution of (<ref>), then $x^*$ will be an optimal solution of (<ref>) and $H(x^*) = F(x^*, \zeta ^*, \beta ^*)$ [28].

As a result, the optimization problem in section <ref>(A) can be written as the following problem:

\begin{equation}
\text{minimize } f(P, \zeta, \beta) = \frac{1}{2} [\sum_{i=1}^{m} (\zeta_i (P_{i,g}^t)^2 + \beta_i \eta_{i,g}^2)] \Delta t
\label{eq:System_Efficiency:Eq19}
\end{equation}

subject to

\begin{equation}
\zeta_i \beta_i \ge 1,
\label{eq:System_Efficiency:Eq20}
\end{equation}

\begin{equation}
(\zeta, \beta) > 0,
\label{eq:System_Efficiency:Eq21}
\end{equation}

§.§ Convex Relaxation Technique

The convex optimization methods can only be applied to problems where the objective function and all constraints are finite and convex. Although the reformulated problem in section <ref>(A) has a convex objective function for a fixed set of $\eta_i$ and $\beta_i$, several constraints still have nonlinearity. In this section, the nonlinear power flow constraints (<ref>), (<ref>), and line flow constraints (<ref>),(<ref>) are replaced with the following linear and quadratic constraints:

\begin{equation}
P_{i,\mathrm{inj}} = \sqrt{2}u_i G_{ii} + \sum_{k \in B} \left(G_{ik}W_{R_{ik}} + B_{ik}W_{I_{ik}}\right),
\label{eq:System_Efficiency:Eq22}
\end{equation}

\begin{equation}
Q_{i,\mathrm{inj}} = -\sqrt{2}u_i B_{ii} + \sum_{k \in B} \left(G_{ik}W_{I_{ik}} - B_{ik}W_{R_{ik}}\right),
\label{eq:System_Efficiency:Eq23}
\end{equation}

\begin{equation}
P_{ik} = \sqrt{2}u_i G_{ik} - \left(G_{ik}W_{R_{ik}} + B_{ik}W_{I_{ik}}\right),
\label{eq:System_Efficiency:Eq24}
\end{equation}

\begin{equation}
Q_{ik} = -\sqrt{2}u_i B_{ik} + \left(B_{ik}W_{R_{ik}} - G_{ik}W_{I_{ik}}\right).
\label{eq:System_Efficiency:Eq25}
\end{equation}

\begin{equation}
W_{R_{ik}}^2 + W_{I_{ik}}^2 \leq 2u_i u_k.
\label{eq:System_Efficiency:Eq26}
\end{equation}

\begin{equation}
\theta_{ik} = \tan^{-1} \left( \frac{W_{I_{ik}}}{W_{R_{ik}}} \right).
\label{eq:System_Efficiency:Eq27}
\end{equation}

Here, (<ref>), (<ref>) are the linear real and reactive power flow equations and (<ref>), (<ref>) are linear line flow equations. The relationship between the convex variables $u_i, W_{I_{ik}}$, and $W_{R_{ik}}$ are defined by the equations (<ref>) and (<ref>). Since (<ref>) is still nonlinear, Taylor series expansion can be used to linearize the equation:

\begin{equation}
\begin{split}
\tan^{-1}\frac{W_{I_{ik}}^{(q)}}{W_{R_{ik}}^{(q)}}=\theta_{ik}+\frac{W_{I_{ik}}^{(q)}}{W_{R_{ik}}^{(q)2}+W_{I_{ik}}^{(q)2}} W_{R_{ik}}- \\ \frac{W_{R_{ik}}^{(q)}}{W_{R_{ik}}^{(q)2}+W_{I_{ik}}^{(q)2}}W_{i_{ik}}.
\label{eq:System_Efficiency:Eq28}
\end{split}
\end{equation}

where, the higher order terms are neglected, and $(W_{I_{ik}}^{(q)}, W_{R_{ik}}^{(q)})$ are the initial estimation. A detailed description of this technique can be found in [31].

The final problem formulation will be as follows:

minimize $f(x, \zeta, \beta)$

subject to

(<ref>), (<ref>), (<ref>)-(<ref>), (<ref>), (<ref>), (<ref>) and (<ref>).

where, $P_{i,inj}, Q_{i,inj}, P_{ik}$, and $Q_{ik}$ are defined by (<ref>), (<ref>), (<ref>), and (<ref>), respectively.

§ SOLUTION TECHNIQUE

In this section, an iterative method is described to solve the reformulated problem.
For a fixed set of $(\zeta, \beta)$, let us consider the following sub-problem of (<ref>):

\begin{equation}
\text{minimize } F(x; \zeta, \beta) = \frac{1}{2} \sum_{i=1}^{m} [\zeta_i (A_i(x))^2 + \beta_i (B_i(x))^2]
\label{eq:System_Efficiency:Eq29}
\end{equation}

Equation (<ref>) can be solved using any standard convex optimization technique. If the optimal solution for (<ref>) is $x^*(\zeta, \beta)$ , then for a fixed set of $x^*$ (<ref>) is reduced to the following problem of 2m variables $(\zeta, \beta)$:

\begin{equation}
\text{minimize } F_{aux}(\zeta, \beta) = \frac{1}{2} \sum_{i=1}^{m} [\zeta_i (A_i(x^*))^2 + \beta_i (B_i(x^*))^2]
\label{eq:System_Efficiency:Eq30}
\end{equation}

subject to,

$\zeta_i \beta_i \ge 1$
$(\zeta, \beta) > 0$

where, $F_{aux}$ denotes the auxiliary problem of $F$. Equation (<ref>) and (<ref>) are solved iteratively until convergence is achieved.

The following algorithm is used to solve the fractional optimization problem:

Fractional Programming Solution Steps

$ \textbf{Step 0: }$ Set i = 0.
$ \textbf{Step 1: }$ Find an optimal solution, $x^*$ for (<ref>) for a fixed $(\zeta, \beta)$.
$ \textbf{Step 2: }$ Find an optimal solution $(\zeta_i, \beta_i)$ for (<ref>) using convex minimization method.
$ \textbf{Step 4: }$ Update the feasible region of (<ref>) by including cutting function $l_i$, where $l_i = 2-\beta_i \sqrt\frac{\zeta_i}{\beta_i} - \zeta_i\sqrt\frac{\beta_i}{\zeta_i} \le0$
$ \textbf{Step 5: }$ Let $i=i+1$ and return to step 0. Continue until converges.

In this paper, the MOSEK optimization toolbox is used to solve the formulated problem [32]. The optimization problem is transformed into the conic quadratic format to fit into MOSEK. MOSEK supports two types of quadratic cones:

* General quadratic cones:

$Q^n = [x \in R^n: x_0\ge\sqrt{\sum_{j=1}^{n-1} x_j^2}]$

* Rotated Quadratic cones:

$Q_r^n = [x \in R^n: 2x_0 x_1\ge\sum_{j=2}^{n-1} x_j^2], x_0 \ge 0, x_1 \ge 0$

All the quadratic parts of the reformulated minimization problem in section <ref>(C) are replaced with rotated quadratic cones and corresponding linear equations. As a result, the transformed problem can be solved using the MOSEK solver.

§ CASE STUDIES

In this literature, the proposed system efficiency model is tested with a notional 12-bus MVAC ship system, a modified IEEE 30-bus, and an IEEE 118-bus system. Each system consists of multiple generators, distribution lines, ESS, and loads at different buses. The load data is generated using the demand pattern of Real-time Dashboard, NYISO [33]. The load profile for different systems can be observed in Fig.<ref>. This paper considers a 24-hour time horizon while solving the optimization problem, where each time step is 1 hour. However, any time horizon and length of time step can be selected based on the system requirement.

24-hours load profile of different systems

§.§ Notional 12-Bus MVAC Ship System

The notional four-zone 12 bus MVAC ship system [34] (shown in Fig. <ref>) consists of two main gas turbine generators (MTG) and two auxiliary gas turbine generators (ATG). The generator parameters can be observed from TABLE <ref>. The system also has 8 ESS (each zone has 2 ESS) and multiple loads, including 2 propulsion motor modules (PMM) at buses 6 and 7 and AC load centers (ACLC) at buses 1,2,3,4,9,10,11, and 12. The ESS data are shown in TABLE <ref>.

Generator Data for Notional MVAC Ship System

Capacity (MW)
Efficiency curve coefficient
Number of Units

ATG     4.7     $a_i = -.133, b_i = .311, c_i = .174$     2

MTG     35     $a_i = -.133, b_i = .311, c_i = .204$     2

The proposed optimization model was run for the notional MVAC ship system; the simulation time was 34.414s, taking 178 iterations to converge. The output power generation can be observed in Fig. <ref>. The system is likely to dispatch the ATGs more than the MTGs to improve the overall efficiency since the capacities of the ATGs are lower than the MTGs.
The ESS charging and discharging schedule and the SOC of the ESS are shown in Figs. <ref> and <ref>, respectively.

A notional MVAC ship system model

Generator output for notional MVAC ship system: (a) MPGM 1, (b) MPGM 2, (c) APGM 1, (d) APGM 2.

ESS charging and discharging schedule for notional MVAC ship system; (a) Zone 1, (b) Zone 2, (c) Zone 3, (d) Zone 4.

SOC of ESS for notional MVAC ship system

ESS data for notional MVAC ship system, IEEE 30-bus, and IEEE 118-bus system

Number of ESS
Capacity (MWh)
Maximum Charging/discharging rate (MW/h)
Miniumum SOC (%)

MVAC ship system     8     2.2     10     20

IEEE 30-bus system     6     15     5     20

IEEE 118-bus system     16     15     5     20

§.§ IEEE 30-bus System

The IEEE 30-bus system [35] has 6 generator buses, 16 load buses, and 42 transmission lines. In addition, the system is modified by including six energy storage systems (ESS) at buses 5, 11, 15, 19, 23 and 27.

The proposed model was then tested for the IEEE 30-bus system, and the observed convergence time was 91.728s for 476 iterations. The generator output schedule is shown in Fig. <ref>. Although the system has six generators, only three produced output during the simulation. Generator graphs with zero output are not included in the figure.
The state of charge (SOC) of the ESS is shown in Fig. <ref>.

Generator output for IEEE 30-bus system at (a) bus 5, (b) bus 11, (c) bus 13.

State of charge (soc) for IEEE 30-bus system

§.§ IEEE 118-bus System

The IEEE 118-bus system [36] has 21 generator buses, 113 load buses, and 179 transmission lines. In addition, the system is modified by including 16 energy storage systems (ESS).

In this study, the IEEE 118-bus system is the largest system where the system efficiency model is applied. The model was run successfully with a reasonable convergence time of 316.22s. Fig. <ref> represents the output power of the generators. For the IEEE 118-bus system, only 9 generators generated power. Only the generators with output are shown in the figure.

Generator output for IEEE 118-bus system

State of charge (SOC) for IEEE 118-bus system

Convergence curves of both objective and tolerance values are shown in  Fig. <ref>. The number of iterations required to converge depends on the number of variables in the system. Since the IEEE 118-bus system has the highest number of variables among the tested systems, it takes the most iteration to converge.

Convergence curves of (a) Objective value (b) Tolerance value

The summary of the results for all systems is shown in TABLE <ref>. The simulations were run on Intel Core i7-10700 CPU, 2.90 GHz processor with 32.0 GB RAM.

Summary results for notional MVAC Ship System, IEEE 30-bus, and IEEE 118-bus system

Number of Iteration
Simulation Time (s)
Fuel Consumption (L)

MVAC ship system     172     34.414     $5.91\times10^4$

IEEE 30-bus system     240     91.728     $3.53\times10^5$

IEEE 118-bus system     581     316.221     $2.28\times10^6$

§.§ Performance Comparison

In this section, a comparative analysis is conducted of the performance of the proposed fuel efficiency model in two distinct domains:

* Comparison in terms of convergence time: This comparison indicates that the model proposed in this paper can be solved efficiently within a reasonable convergence time. As a result, the model can be applied to large-scale systems where the nonlinear programming model takes excessive time to converge.

* Comparison in terms of fuel consumption: The proposed model consumes significantly less amount of fuel compared to the other models with different generator dispatches. The results indicate that the proposed model is the most efficient and optimal.

§.§.§ Comparison in terms of convergence time

The convergence time of NLP and FP models is compared in this subsection. The nonlinear optimization model from section <ref>(A) is solved with MATLAB NLP (fmincon function) for the notional MVAC ship system. The solution took an extensive time (more than 8 hours) to converge for the notional MVAC ship system with 24-time steps where the convergence time of the FP model was only 34.414s for the same system. As a result, the NLP makes the solution procedure impractical to apply to extensive systems. Moreover, the fuel consumed during the operation was $6.53\times10^4$L, higher than the fuel consumed with the proposed FP model. The FP model is clearly more advantageous than the nonlinear programming optimization, even with a higher number of variables.
The performance comparison between FP and NLP is shown in TABLE <ref>.

Performance comparison between FP and NLP for notional MVAC ship system

System     Optimization Model     Total Variables     Simulation Time (s)     Fuel Consumption (L)

25emNotional MVAC Ship System     FP     2592     34.414     $5.91\times10^4$

NLP     1488     30233.884     $6.53\times 10^4$

2*IEEE 30-bus     FP     6960     91.728     $3.53\times10^5$

NLP     4128     did not converge     NA

2*IEEE 118-bus     FP     28128     316.221     $2.28\times10^6$

NLP     17088     did not converge     NA

§.§.§ Comparison in terms of fuel consumption

The proposed system efficiency model is compared with three other models (as listed in TABLE <ref>) to demonstrate the fuel efficiency. The difference between the models is in the generator dispatch allowed during the simulation. The number and type of generators allowed for each model during the operation are indicated by the 'Dispatch' column.

The comparison results are shown in Fig. <ref>. It can be observed that the proposed model has the lowest fuel consumption among all models. Initially, all models consume almost similar amounts of fuel. However, the difference in fuel consumption increases with the number of time steps. This observation highlights the superior efficiency of the proposed model in comparison to the models examined in this section.

Fuel consumption comparison between different models of notional MVAC system

Summary of generator dispatch of different models of notional MVAC ship system

Proposed Model     2 MTGs and 2 ATGs

Model A     2 MTGs

Model B     1 MTG and 1 ATG

Model C     Equal power-sharing

§ CONCLUSION

This study has addressed the challenge of the fuel consumption minimization problem to enhance the system efficiency and reduce the operating cost of the power generation units. The traditional approaches typically focus on maximizing the efficiency function or minimizing the generator cost function to achieve optimal fuel consumption for the system. However, these approaches do not account for the fuel consumption rate directly and are impractical to implement in real-world systems where optimizing fuel use is the objective. In addition, existing studies that have used the fuel consumption curve to formulate the optimization problem have numerous limitations, including incompatibility to apply to large AC systems.

As a result, it is crucial to incorporate a function that directly represents fuel consumption to enhance the system's fuel efficiency. This study introduced a novel objective function based on a sum-of-ratios approach, providing a straightforward representation of the fuel consumption rate. The sum-of-ratios problem was effectively solved by leveraging a fractional programming (FP) reformulation technique, resulting in successful fuel consumption minimization. Moreover, the low convergence time of the solution makes the model suitable for large-scale systems. While the model stands out in its uniqueness and effectiveness compared to other approaches, future research will concentrate on implementing a distributed algorithm to enhance scalability for larger and more complex systems.

§ ACKNOWLEDGEMENT

The information, data, or work presented herein was partly funded by the U.S. Office of Naval Research under the award numbers N000142212239 and N000142112124.

[1]
K. D. Mistry and R. Roy, “Enhancement of loading capacity of distribution
system through distributed generator placement considering techno-economic
benefits with load growth,” vol. 54, pp. 505–515. [Online]. Available:

[2]
Y. Yao and Q. Gao, “Efficiency optimization of a power station with different
generators,” vol. 145, pp. 440–445. [Online]. Available:

[3]
S. Orero and M. Irving, “A genetic algorithm modelling framework and solution
technique for short term optimal hydrothermal scheduling,” vol. 13, no. 2,
pp. 501–518, conference Name: IEEE Transactions on Power Systems.

[4]
X. Yuan, Y. Zhang, L. Wang, and Y. Yuan, “An enhanced differential evolution
algorithm for daily optimal hydro generation scheduling,” vol. 55, no. 11,
pp. 2458–2468. [Online]. Available:

[5]
T. V. Vu, D. Perkins, D. Gonsoulin, C. S. Edrington, B. Papari, K. Schoder,
M. Stanovich, and M. Steurer, “Large-scale distributed control for MVDC
ship power systems,” in IECON 2018 - 44th Annual Conference of the
IEEE Industrial Electronics Society.1em plus 0.5em minus
0.4emIEEE, pp. 3431–3436. [Online]. Available:

[6]
C. Hernandez-Aramburo, T. Green, and N. Mugniot, “Fuel consumption
minimization of a microgrid,” vol. 41, no. 3, pp. 673–681, conference Name:
IEEE Transactions on Industry Applications.

[7]
M. Lundh, W. Garcia-Gabin, K. Tervo, and R. Lindkvist, “Estimation and
optimization of vessel fuel consumption,” vol. 49, no. 23, pp. 394–399.
[Online]. Available:

[8]
F. D. Kanellos, G. J. Tsekouras, and N. D. Hatziargyriou, “Optimal demand-side
management and power generation scheduling in an all-electric ship,” vol. 5,
no. 4, pp. 1166–1175, conference Name: IEEE Transactions on Sustainable

[9]
S. Barsali, C. Miulli, and A. Possenti, “A control strategy to minimize fuel
consumption of series hybrid electric vehicles,” vol. 19, no. 1, pp.
187–195, conference Name: IEEE Transactions on Energy Conversion.

[10]
A. Accetta and M. Pucci, “Energy management system in DC micro-grids of
smart ships: Main gen-set fuel consumption minimization and fault
compensation,” vol. 55, no. 3, pp. 3097–3113, conference Name: IEEE
Transactions on Industry Applications.

[11]
K. Satpathi, V. M. Balijepalli, and A. Ukil, “Modeling and real-time
scheduling of DC platform supply vessel for fuel efficient operation,”
vol. 3, no. 3, pp. 762–778, conference Name: IEEE Transactions on
Transportation Electrification.

[12]
B. Chowdhury and S. Rahman, “A review of recent advances in economic
dispatch,” vol. 5, no. 4, pp. 1248–1259, conference Name: IEEE
Transactions on Power Systems.

[13]
A. B. Kunya, A. S. Abubakar, and S. S. Yusuf, “Review of economic dispatch in
multi-area power system: State-of-the-art and future prospective,” vol. 217,
p. 109089. [Online]. Available:

[14]
R. W. Freund and F. Jarre, “Solving the sum-of-ratios problem by an
interior-point method.”

[15]
T. Kuno, “A branch-and-bound algorithm for maximizing the sum of several
linear ratios.”

[16]
J. G. Carlsson and J. Shi, “A linear relaxation algorithm for solving the
sum-of-linear-ratios problem with lower dimension,” vol. 41, no. 4, pp.
381–389. [Online]. Available:

[17]
A. Charnes and W. W. Cooper, “Programming with linear fractional
functionals,” vol. 9, no. 3, pp. 181–186, _eprint:
https://onlinelibrary.wiley.com/doi/pdf/10.1002/nav.3800090303. [Online].

[18]
S. Schaible, “Parameter-free convex equivalent and dual programs of fractional
programming problems,” vol. 18, no. 5, pp. 187–196. [Online]. Available:

[19]
W. Dinkelbach, “On nonlinear fractional programming,” vol. 13, no. 7, pp.
492–498. [Online]. Available:

[20]
M. Toloo, “An equivalent linear programming form of general linear fractional
programming: A duality approach,” vol. 9, no. 14, p. 1586. [Online].
Available: <https://www.mdpi.com/2227-7390/9/14/1586>

[21]
M. Jayalakshmi, “A new approach for solving quadratic fractional programming

[22]
Department of Mathematical Engineering, Yildiz Technical University, Istanbul,
Turkey, M. Sivri, I. Albayrak, Department of Mathematical Engineering,
Yildiz Technical University, Istanbul, Turkey, G. Temelcan, and Department
of Computer Programming, Istanbul Aydin University, Istanbul, Turkey, “A
novel approach for solving quadratic fractional programming problems,”
vol. 9, no. 2, pp. 199–209. [Online]. Available:

[23]
Y. Almogy and O. Levin, “A class of fractional programming problems,”
vol. 19, no. 1, pp. 57–67, publisher: INFORMS. [Online]. Available:

[24]
J. E. Falk and S. W. Palocsay, Optimizing the Sum of Linear Fractional
Functions.1em plus 0.5em minus 0.4emPrinceton University
Press, pp. 221–258. [Online]. Available:

[25]
Y. Jong, “Practical global optimization algorithm for the sum-of-ratios

[26]
J.-F. Tsai, “Global optimization of nonlinear fractional programming problems
in engineering design,” vol. 37, no. 4, pp. 399–409. [Online]. Available:

[27]
K. Shen and W. Yu, “Fractional programming for communication systems—part i:
Power control and beamforming,” vol. 66, no. 10, pp. 2616–2630. [Online].
Available: <https://ieeexplore.ieee.org/document/8314727/>

[28]
H. Konno, T. Kuno, and Y. Yajima, “Global minimization of a generalized convex
multiplicative function,” vol. 4, no. 1, pp. 47–62. [Online]. Available:

[29]
S. Schaible and J. Shi, “Fractional programming: The sum-of-ratios case,”
vol. 18, no. 2, pp. 219–229. [Online]. Available:

[30]
Battery minimum state of charge. [Online]. Available:

[31]
M. I. Anam, T. Vu, K. Murari, C. Schegan, and M. Steurer, “System efficiency
model for microgrids using fractional programming,” in 2023 IEEE
Electric Ship Technologies Symposium (ESTS).1em plus 0.5em minus
0.4emIEEE, pp. 36–41. [Online]. Available:

[32]
Mosek - downloads. [Online]. Available: <https://www.mosek.com/downloads/>

[33]
Real-time dashboard. [Online]. Available:

[34]
ESRDC  electric ship research and development consortium).
[Online]. Available:

[35]
IEEE 118-bus system - illinois center for a smarter electric grid (ICSEG).
[Online]. Available:

[36]
IEEE 30-bus system - illinois center for a smarter electric grid (ICSEG).
[Online]. Available: <https://icseg.iti.illinois.edu/ieee-30-bus-system/>


# Quasinormal modes in semi-open systems

Leonardo Solidoro<EMAIL_ADDRESS>School of Mathematical
Sciences, University of Nottingham, University Park, Nottingham, NG7 2RD, UK
Nottingham Centre of Gravity, University of Nottingham, University Park,
Nottingham NG7 2RD, UK  Sam Patrick  Department of Physics, King’s College
London, University of London, Strand, London, WC2R 2LS, UK  Ruth Gregory
Department of Physics, King’s College London, University of London, Strand,
London, WC2R 2LS, UK  Perimeter Institute, 31 Caroline Street North, Waterloo,
ON, N2L 2Y5, Canada  Silke Weinfurtner  School of Mathematical Sciences,
University of Nottingham, University Park, Nottingham, NG7 2RD, UK  Nottingham
Centre of Gravity, University of Nottingham, University Park, Nottingham NG7
2RD, UK  Centre for the Mathematics and Theoretical Physics of Quantum Non-
Equilibrium Systems, University of Nottingham, Nottingham, NG7 2RD, UK

###### Abstract

Astrophysical black holes are open systems which, when perturbed, radiate
quasi-normal modes (QNMs) to infinity. By contrast, laboratory analogues are
necessarily finite-sized, presenting a potential obstacle to exciting QNMs in
the lab. In this study, we investigate how the QNM spectrum of a toy-model
black hole is modified when the system is enclosed by a partially reflecting
wall. Counter to expectation, we demonstrate that QNMs not only persist in
finite-sized systems, but the number of detectable modes increases.
Furthermore, we show that QNMs in this set-up can be easily excited by
incoherent background noise. Our findings align with studies exploring the
spectral stability of black holes, such as those examining small modifications
of the surrounding gravitational field. Importantly, our work paves the way
for exploring spectral stability in laboratory systems, enabling experimental
investigation of black hole features in closed systems.

_Introduction.—_ The success of the black hole spectroscopy program is founded
on the use of gravitational waves produced by merging binaries as a probe of
the system’s properties [1]. At late times during the so-called _ringdown_
phase, the remnant emits a characteristic signal as it sheds energy and
relaxes into its final state. This signal consists of a discrete set of damped
harmonic oscillators referred to as quasinormal modes (QNMs) [2, 3, 4, 5, 6],
which are natural resonances of linear open systems. Next generation
gravitational wave detectors such as LISA and Einstein telescope will profit
from enhanced sensitivity [7, 8], making it paramount to understand how the
ringdown is influenced by various high-order effects, e.g. environmental
factors [9] and nonlinearities [10].

Recent investigation into the spectral stability of black holes has
demonstrated the susceptibility of the QNM spectrum to small, local
perturbations to the scattering potential experienced by gravitational waves
around a black hole [11, 12, 13, 14]. In some cases, the spectrum can even
admit additional long-lived modes called quasi-bound states, potentially
leading to higher-order effects as a result of their long lifetimes [15, 16].
Understanding such modifications of the QNM spectrum may have important
consequences for the black hole spectroscopy programme [17, 18]. In this
regard, it is natural to consider non-perturbative scenarios where energy is
retained by the system due to an effective reflection mechanism acting on
escaping radiation. Such situations arise naturally in gravity in the context
of neutron stars [4], exotic compact objects [19, 20, 21, 22] and
asymptotically anti-de Sitter spacetimes [6].

Gravity simulators present an ideal platform where such an enquiry can be
pursued. These are physical systems in which excitations behave as though they
propagate on a curved spacetime, permitting the experimental investigation of
elusive phenomena like Hawking radiation [23, 24, 25, 26] and rotational
superradiance [27, 28, 29, 30, 31]. Notably, ringdown signals have been
studied in polariton superfluids [32], optical solitons [33] and in
hydrodynamic systems [34, 35]. Although treating these platforms as
effectively open systems is usually justified by a combination of energy
dissipation, limited temporal evolution and engineering of absorptive
boundaries, a complete analysis should not neglect finite-size effects
resulting from confinement mechanisms. This is particularly true of our recent
superfluid helium gravity simulator [36], where dissipation is strongly
reduced and many reflections from the container wall can occur during the
experimental timeframe.

All of this serves to emphasise the importance of understanding how finite-
size effects modify the QNM spectrum in the context of both black hole physics
and gravity simulators. In this Letter, we study the effect of a partially
reflective boundary on the QNM spectrum of the Pöschl-Teller (PT) potential, a
model which appears in both black hole physics [6, 2] and analogue systems
[37, 33], demonstrating how finite-size effects endow the system with an
additional control parameter that can be exploited in the experimental study
of QNMs.

_General set-up.—_ A simple system which captures the essential physics of
QNMs is the one-dimensional wave equation with a repulsive potential $V(x)$
[4, 6],

$\displaystyle\left(\frac{\partial^{2}}{\partial
t^{2}}-\frac{\partial^{2}}{\partial x^{2}}+V(x)\right)\Psi=0\,,$ (1)

where the wave speed is set to unity and we assume $V\to 0$ as
$x\to\pm\infty$. QNMs are resonances of open systems, i.e. systems that lose
energy, in this case, to spatial infinity. By looking for solutions of the
form $\Psi(x,t)=e^{-i\omega t}\psi(x)$, the wave equation reduces to a
Schrödinger-like equation
$\psi^{\prime\prime}+\left[\omega^{2}-V(x)\right]\psi=0\,,$ and the set of
QNMs are found as the set of frequencies whose related modes satisfy the
purely out-going boundary conditions,

$\displaystyle\lim_{x\to\pm\infty}\psi\propto e^{\pm i\omega x}.$ (2)

More generally, the problem of calculating QNMs can be formulated as a
scattering problem. Consider a solution of (1) consisting of an incident and a
reflected contribution on the right and a transmitted one on the left,

$\displaystyle\psi(x)=\begin{cases}A_{i}\,e^{-i\omega x}+A_{r}\,e^{+i\omega
x}\quad&,\,\,x\to+\infty\,,\\\ A_{t}\,e^{-i\omega
x}\quad&,\,\,x\to-\infty\end{cases}\,.$ (3)

From this, one can define the reflection and transmission coefficients
$R(\omega)=A_{r}/A_{i}$ and $T(\omega)=A_{t}/A_{i}$ respectively, and one sees
that the QNM frequencies determined by (2) are poles of the reflection
coefficient [6].

Consider now the semi-open problem where a wall of reflectivity $\varepsilon$
is situated at $x_{b}$ such that $V(x_{b})$ is effectively zero. The boundary
condition becomes,

$\displaystyle\psi(x\to x_{b})\propto\varepsilon
e^{i\omega(x-x_{b})}+e^{-i\omega(x-x_{b})}\,.$ (4)

We choose a constant $\varepsilon\in[0,1]$, allowing us to vary the amount of
absorption at the wall: for $\varepsilon=1$ the wall is fully reflective
whilst for $\varepsilon=0$ we recover the open problem 111Note that our choice
of $\varepsilon$ real and positive corresponds to a Neumann boundary condition
in the case of a complete reflection, whereas choosing $\varepsilon$
negative/complex would lead to a Dirichlet/Robin boundary condition. In
general, $\varepsilon$ could depend on $\omega$ but for the sake of
simplicity, we model it here as a constant.. The compatibility of (4) with (3)
implies the following resonance condition [39],

$\displaystyle\varepsilon R(\omega)e^{2i\omega x_{b}}=1\,.$ (5)

The QNM condition is retrieved in the $\varepsilon\to 0$ limit, whilst for a
finite reflectivity one should expect the resonances to deviate from the ones
of the open system. This is well understood in the context of spectral
stability analysis, where a small perturbation of $V$ can significantly modify
the spectrum [17]. The continuous evolution of the resonances in the complex
plane as $\varepsilon$ increases can be described in terms of integral curves
of an ordinary differential equation [40]. The resonance condition (5) can be
expressed as $\varepsilon-R^{-1}(\omega_{n}^{\varepsilon})=0$, where
$\omega_{n}^{\varepsilon}$ are the resonant frequencies indexed by
$n=0,1,...$, and the superscript indicates the value at a particular
$\varepsilon$. By differentiating with respect to the reflectivity, one gets,

$\displaystyle\frac{\partial\omega_{n}}{\partial\varepsilon}=-\frac{R^{2}}{R^{\prime}}\bigg{|}_{\omega=\omega_{n}}\equiv\mathrm{H}(\omega_{n})\,,$
(6)

where $R^{\prime}=\partial_{\omega}R$. The right-hand side defines a vector
field along which the real and imaginary parts of the resonances evolve as we
change the boundary’s reflectivity.

_Pöschl-Teller potential.—_ We now apply this framework to the Pöschl-Teller
(PT) potential,

$\displaystyle V_{\mathrm{pt}}(x)=\frac{V_{0}}{\cosh^{2}(\alpha\,x)}\,.$ (7)

In this case, (1) has an analytical solutions in the form $e^{-i\omega
t}\psi(x)$. Defining a new variable $y=1/(1+e^{2\alpha x})$, the spatial part
can be written as,

$\displaystyle\psi(y)$
$\displaystyle=(1-y)^{-i\omega/2\alpha}\left[A\,y^{-i\omega/2\alpha}{}_{2}F_{1}(a,b,c;y)+\right.$
$\displaystyle\left.+B\,y^{i\omega/2\alpha}{}_{2}F_{1}(1+a-c,1+b-c,2-c;y)\right]\,,$
(8)

where $A,B$ are constants and ${}_{2}F_{1}$ hypergeometric functions [41]
whose parameters are given by,

$a,b=\frac{1}{2}-i\omega/\alpha\pm\frac{1}{2}\sqrt{1-\frac{4V_{0}}{\alpha^{2}}},\qquad
c=1-i\omega/\alpha\,,$ (9)

with the $+$ ($-$) sign for $a$ ($b$). Without loss of generality, from now on
we take $V_{0}=1$. To impose the resonance condition (5), we need to define
the reflection coefficient. Assuming that $V(x_{b})$ is close to vanishing,
(Quasinormal modes in semi-open systems) is well-approximated by its
asymptotic expansion and the in- and out-going wave contributions can be
identified (see Appendix B). Their ratio gives the reflection coefficient,

$\displaystyle
R=\frac{\Gamma(a)\Gamma(b)}{\Gamma(c-a)\Gamma(c-b)}\frac{\Gamma(i\omega/\alpha)}{\Gamma(-i\omega/\alpha)}.$
(10)

Poles of this expression give the QNMs of the open system,

$\displaystyle\omega^{\varepsilon=0}_{n}$
$\displaystyle=\pm\sqrt{1-\alpha^{2}/4}-i\alpha(2n+1)/2\,,$ (11)

where $n$ is called the overtone number. All QNM frequencies of the PT
potential share the same real part and differ only in their imaginary part.
The mode with $n=0$ is called the fundamental mode, since it has the longest
lifetime, whereas overtones with $n>0$ are characterised by stronger damping
[6]. This particular feature, also relevant for astrophysical black holes,
renders the observation of overtones difficult. Consequently, the analysis of
the systems’ properties is often limited to the detection of the fundamental
mode.

Figure 1: a. Trajectories in the complex plane of the first 11 resonances of
the PT potential with a partially-reflective barrier relative to the
fundamental QNM
$\omega_{n=0}^{\varepsilon=0}\equiv\tilde{\omega}_{r}+i\tilde{\omega}_{i}$.
Here we take $\alpha=0.4$ and the position of the wall is fixed at $x_{b}=10$.
The large yellow dots represent QNMs of the open system, whereas the blue
lines are obtained varying $\varepsilon$. Coloured markers label the
resonances at particular values of reflectivity. b. Trajectories of the
resonances as reflectivity increases, for $\alpha=0.5$ (Green, dashed), $=0.4$
(Blue, solid), $=0.3$ (Orange, dotted) and $x_{b}=10$. The yellow dots are
$\omega_{n=0,1}^{\varepsilon=0}$ and the lines end at $\varepsilon=1$. c. The
migration of the resonance is compared against the vector field
$\mathrm{H}(\omega)$ for $\alpha=0.4$. The red cross is placed on the
repelling point, i.e. the pole of $\mathrm{H}$.

Resonances of the semi-open system are obtained by numerically solving (5). In
Fig. 1a, we display their continuous migration in the complex plane as
$\varepsilon$ is varied. Large yellow dots represent QNMs of the open system
for $\varepsilon=0$ and, as $\varepsilon$ is increased, the resonances migrate
down toward the real line. We can identify three qualitatively different types
of QNM when $\varepsilon=1$. Modes with $\mathrm{Re}[\omega]=0$ correspond to
purely damped solutions (discussed in [40]). Modes with
$0<\mathrm{Re}[\omega]<\tilde{\omega}_{r}$ are damped oscillatory modes which
are localised between the hard wall and the barrier $V$. These are quasibound
states (QBS) which need to tunnel through $V$ to escape to $x\to-\infty$,
leading to long lifetimes. Finally, modes with
$\mathrm{Re}[\omega]>\tilde{\omega}_{r}$ pass over the barrier and experience
only a small amount of reflection due to the inhomogeneous nature of $V$.
Consequently, they are more strongly damped and can be identified with the new
overtones of the system. These overtones possess distinctly different real
frequencies and an imaginary part which is reduced significantly relative to
the open system.

In Fig. 1b, we provide a closeup of the $n=0,1$ modes and show the effect of
changing $\alpha$, i.e. the width of the PT potential. As $\alpha$ gets
larger, the potential barrier becomes narrower and more QBSs can fit in the
cavity formed by the hard wall and potential barrier. In the complex $\omega$
plane, one would then see all the $\omega^{\varepsilon=1}_{n}$ moving toward
lower $\mathrm{Re}[\omega]$ and the trajectories emerging from
$\omega^{\varepsilon=0}_{n}$ orbiting around the yellow points in Fig. 1a in
the clockwise direction. However, at critical values of $\alpha$, there can be
discontinuous jumps where the migrations lines ending at
$\omega^{\varepsilon=1}_{n}$ change which of the $\omega^{\varepsilon=0}_{n}$
they are connected to. These discontinuous jumps are related to repelling
points of (5), illustrated as the red cross in Fig. 1c. A migration line is an
integral curve of the vector field $\mathrm{H}(\omega)$, defined in (6), and
is deflected either to the left or to the right of the repelling point, i.e.
the poles of $\mathrm{H}(\omega)$, as the resonances move down toward the real
axis. A small change in $\alpha$ can change the vector field such that a left-
directed line suddenly changes to a right-directed line. This is exemplified
by the transition from $\alpha=0.4$ to $\alpha=0.5$ on Fig. 1b where for
$\alpha=0.4$, the line emerging from $\omega^{\varepsilon=0}_{n=0}$
($\omega^{\varepsilon=0}_{n=1}$) is directed left (right), but switches
direction when increasing to $\alpha=0.5$. Although the trajectories are
discontinuous when varying $\alpha$, and similar behaviour would occur when
varying $x_{b}$, what matters for us is that for $\alpha$ and $x_{b}$ fixed,
varying $\varepsilon$ defines a continuous map from the open system resonances
to the QBS and the new overtones of the semi-open system. Hence, $\varepsilon$
is a dial which allows us to smoothly switch on this additional structure.

Figure 2: a. Ringdown signal in time and frequency domain for the same $x_{b}$
and $\alpha$ in Fig. 1a. The logarithm of the absolute value of the signal
close to the reflective wall is shown on the left panels. For the open-system
$\varepsilon=0$, we observe the characteristic QNM response whilst for
$\varepsilon\neq 0$, the signal is characterised by a fast (slowly) decaying
signal at early (late) times. b. We show Fourier transform the of the signal
at $\bar{x}_{1}$ near the wall (dotted line) and on the open side of the
system at $\bar{x}_{2}$ (solid). The first overtone (the peak to the right of
$\mathrm{Re}[\omega]>1$) becomes more pronounced for larger $\varepsilon$
since the damping is reduced relative to $\varepsilon=0$, and corresponds to
the short-lived part of the time domain signal. c. We plot the Fourier
transform of the signal at each point to illustrate the spatial features for
the resonant modes against the PT potential (thick white line). The signal is
normalised such that the colour bar goes from 0 (black) to 1 (white). The
dotted and the solid lines label $\bar{x}_{1}$ and $\bar{x}_{2}$ respectively.
The spectrum averaged over the $x$ axis is showed on the right side.

_Numerical simulations.—_ A consequence of the reflective boundary is that the
new overtones have distinct real frequencies and longer lifetimes,
facilitating their detection. We demonstrate this via a numerical simulation.
The wave equation (1) is simulated using a method of lines algorithm (details
in Appendix A). In Fig. 2, a Gaussian pulse is evolved toward the barrier from
the side of the open boundary. We illustrate the response signal, in both the
time and frequency domains (left and right panels respectively) after the
initial pulse has passed for various $\varepsilon$, i.e. as the boundary is
made more reflective, focusing on values which make the decay times of $n=0,1$
comparable (see Fig. 1). For the open system, we observe the characteristic
QNM waveform, corresponding to a singly peaked frequency spectrum. However,
for semi-open systems, the signal becomes more complicated, characterised by
interference between multiple frequencies in the time domain which correspond
to more than one peak in Fourier space. Notably, the peak corresponding to
$\omega^{\varepsilon=0}_{n=0}$ splits in two, the left peak corresponding to
the highest frequency QBS and the right peak to the first of the new
overtones. For a partially absorbing wall, we observe that this overtone is
actually the dominant signal when we measure near the open boundary. When the
wall is fully reflective, the dominant signal comes from the QBS due to their
long lifetimes, although our simulations suggest that the first overtone is
still detectable. In Fig. 2c, we illustrate the full spatial waveform as a
function of frequency, illustrating how the different resonant modes
(indicated by bright horizontal lines) are bound by the potential barrier. As
expected, modes scattering with the wider part of the barrier experience less
leakage. Note that although format in Fig. 2a will be familiar to
gravitational community, the equivalent representations in panels b and c
emerge naturally in the data analysis of analogue experiments [35, 42] and
carry more information about the spatial profile of the modes.

Although Gaussian pulses are frequently used in numerical studies to stimulate
the system [6], a more natural perturbing mechanism from the experimental
perspective is broadband mechanical noise [35, 36]. In Fig. 3, we take
$\Psi(x,t=0)=0$ as initial condition but generate random noise on the side of
the open boundary at each time step for the duration of the simulation. We
then Fourier transform the signal after the noise has propagated across the
system and take the average on the side of the barrier near the wall. We
collect statistics over 500 iterations and display the mean and the variance.
In the open system, there is no distinct signal; at high frequency we observe
the noise level and at lower frequencies a suppression due to the presence of
the barrier. For partial absorption, we observe a series of peaks emerging
above the noise level, which become very pronounced for the fully reflective
wall. Notably, the first overtone appears as a distinct signal in the data,
suggesting that the combination of a finite size system with mechanical noise
may be sufficient ingredients to detect QNM overtones in analogue experiments.
This method of resonant stimulation is akin to that of the driven, damped
harmonic oscillator, where the energy input from the driving force is balanced
by the energy loss due to dissipation (in our case resulting from the open
boundary).

_Discussion.—_ Finite-size constraints of analogue experiments might initially
appear as an obstacle to simulating black hole ringdown - a phenomenon
typically associated with open (i.e. spatially unbounded) systems. Counter to
this expectation, our results demonstrate that the introduction of such
confining mechanisms, in fact, becomes a valuable tool for studying the
quasinormal modes (QNMs) and the role of the overtones in shaping the ringdown
signal. The set-up we have presented (based on a Poschl-Teller potential)
captures the essential physics of experimentally realisable black hole
analogues, where the partially reflective boundary models the experimental
confinement and the open boundary represents the analogue black hole horizon.
As the reflectivity of the closed boundary is increased from zero, the QNMs of
the open system migrate continuously along trajectories in the complex
frequency plane toward the real axis. In addition to the familiar quasibound
states of semi-open systems, the modified spectrum contains a hierarchy of
overtones with distinct oscillation frequencies and substantially longer
lifetimes compared to the open system.

These two features significantly increase the prospect of a detecting QNM
overtones in an analogue experiment, as exemplified by our numerical
simulations. An added consequence is that the frequency response of the semi-
open system to a random driving signal is sharply peaked around the QNMs, a
feature which is absent in the open system as it relies on constructive
interference with reflected waves. A measurement of such a signal could prove
particularly natural for experiments where noise (e.g. thermal or mechanical)
is always present. Indeed, this mechanism may be responsible for the QNM
oscillations measured in [37].

Figure 3: The systems response in the frequency domain when stimulated with
mechanical noise. We take the average over 500 runs (black lines) and display
the standard deviation as the shaded regions. The semi-open system is
characterised by the presence of well defined resonance peaks. In particular,
we manage to resolve the first overtone above the noise level. The data is for
a potential with $\alpha=1$ and a wall at $x_{b}=10$.

From a gravity perspective, the semi-open problem has been addressed in the
context of neutron stars/compact objects and anti-de Sitter spacetimes. The
focus (in the case of the former) is usually on the new bound states that
arise in the spectrum, and (in the latter) on echoes of the usual QNM
response. In our case, where the reflecting boundary is close to the peak of
the scattering potential, we have shown that the new overtones are a prominent
feature of linear response, and that random noise is a viable probe of the
spectrum. In addition to the aforementioned astrophysical settings,
stimulating QNMs in this fashion may also be relevant for a black hole
surrounded by matter [43, 44, 9] or hypothetical ultra-light fields [45, 46].

Furthermore, discerning the roles of nonlinearities and overtones in shaping
the ringdown waveform is a central topic of research in the gravitational wave
community [47, 48]. Understanding the interplay of these two competing effects
is therefore desirable. Finite size black hole simulators may be useful in
this regard, since both the nonlinearities and overtone lifetimes can be
experimentally controlled (the former through signal amplitudes and the latter
using the boundary reflectivity). This feature stands out as an enhanced
capability of semi-open systems compared with gravitational black holes which
radiate to infinity, where a consequence of the no-hair theorem is that the
overtone lifetime is not an adjustable parameter [6].

_Acknowledgements.—_ L.S. and S.P. would like to thank Théo Torres for
providing valuable insight on spectral stability based on unpublished work
with Sam Dolan. We are grateful to Patrik Švančara, Pietro Smaniotto and
Thomas Sotiriou for helpful discussions and comments on the manuscript. We
also thank Paolo Pani for feedback. L.S. and S.W. gratefully acknowledge the
support of the Leverhulme Research Leadership Award (RL-2019-020). S.P., R.G.
and S.W. extend their appreciation to the Science and Technology Facilities
Council for their generous support in Quantum Simulators for Fundamental
Physics (QSimFP), (ST/T005858/1 and ST/T006900/1), as part of the UKRI Quantum
Technologies for Fundamental Physics program. S.W. also acknowledges the Royal
Society University Research Fellowship (UF120112). R.G. acknowledges support
from the Perimeter Institute. Research at Perimeter Institute is supported by
the Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Research, Innovation and Science.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ , Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. 116, 061102 (2016).
  * Ferrari and Mashhoon [1984] V. Ferrari and B. Mashhoon, New approach to the quasinormal modes of a black hole, Physical Review D 30, 295 (1984).
  * Echeverria [1989] F. Echeverria, Gravitational-wave measurements of the mass and angular momentum of a black hole, Phys. Rev. D 40, 3194 (1989).
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Quasi-normal modes of stars and black holes, Living Reviews in Relativity 2, 1 (1999).
  * Cardoso  _et al._ [2009] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Geodesic stability, lyapunov exponents, and quasinormal modes, Phys. Rev. D 79, 064016 (2009).
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Classical and Quantum Gravity 26, 10.1088/0264-9381/26/16/163001 (2009).
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, _et al._ , Laser interferometer space antenna, arXiv preprint arXiv:1702.00786 (2017).
  * Maggiore  _et al._ [2020] M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse, S. Foffa, J. García-Bellido, _et al._ , Science case for the Einstein telescope, J. Cosmo. Astro Phys. 2020, 050 (2020).
  * Cannizzaro  _et al._ [2024] E. Cannizzaro, T. F. Spieksma, V. Cardoso, and T. Ikeda, The impact of plasma on the relaxation of black holes, arXiv preprint arXiv:2405.05315 (2024).
  * Cheung  _et al._ [2023] M. H.-Y. Cheung, V. Baibhav, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, F. Duque, T. Helfer, E. Shukla, and K. W. K. Wong, Nonlinear effects in black hole ringdown, Phys. Rev. Lett. 130, 081401 (2023).
  * Nollert [1996] H.-P. Nollert, About the significance of quasinormal modes of black holes, Phys. Rev. D 53, 4397 (1996).
  * Aguirregabiria and Vishveshwara [1996] J. Aguirregabiria and C. Vishveshwara, Scattering by black holes: a simulated potential approach, Physics Letters A 210, 251 (1996).
  * Nollert and Price [1999] H.-P. Nollert and R. H. Price, Quantifying excitations of quasinormal mode systems, Journal of Mathematical Physics 40, 980 (1999).
  * Cheung  _et al._ [2022] M. H.-Y. Cheung, K. Destounis, R. P. Macedo, E. Berti, and V. Cardoso, Destabilizing the fundamental mode of black holes: The elephant and the flea, Phys. Rev. Lett. 128, 111103 (2022).
  * Witek  _et al._ [2013] H. Witek, V. Cardoso, A. Ishibashi, and U. Sperhake, Superradiant instabilities in astrophysical systems, Phys. Rev. D 87, 043513 (2013).
  * Huang and Zhang [2021] Y. Huang and H. Zhang, Quasibound states of charged dilatonic black holes, Phys. Rev. D 103, 044062 (2021).
  * Jaramillo  _et al._ [2021] J. L. Jaramillo, R. P. Macedo, and L. A. Sheikh, Pseudospectrum and black hole quasinormal mode instability, Phys. Rev. X 11, 031003 (2021).
  * Jaramillo  _et al._ [2022] J. L. Jaramillo, R. P. Macedo, and L. A. Sheikh, Gravitational wave signatures of black hole quasinormal mode instability, Phys. Rev. Lett. 128, 211102 (2022).
  * Maggio  _et al._ [2017] E. Maggio, P. Pani, and V. Ferrari, Exotic compact objects and how to quench their ergoregion instability, Phys. Rev. D 96, 104047 (2017).
  * Testa and Pani [2018] A. Testa and P. Pani, Analytical template for gravitational-wave echoes: Signal characterization and prospects of detection with current and future interferometers, Phys. Rev. D 98, 044018 (2018).
  * Maggio  _et al._ [2019] E. Maggio, V. Cardoso, S. R. Dolan, and P. Pani, Ergoregion instability of exotic compact objects: Electromagnetic and gravitational perturbations and the role of absorption, Phys. Rev. D 99, 064007 (2019).
  * Maggio  _et al._ [2020] E. Maggio, L. Buoninfante, A. Mazumdar, and P. Pani, How does a dark compact object ringdown?, Phys. Rev. D 102, 064053 (2020).
  * Euvé  _et al._ [2016] L.-P. Euvé, F. Michel, R. Parentani, T. G. Philbin, and G. Rousseaux, Observation of noise correlated by the Hawking effect in a water tank, Phys. Rev. Lett. 117, 121301 (2016).
  * Steinhauer [2014] J. Steinhauer, Observation of self-amplifying Hawking radiation in an analogue black-hole laser, Nat. Phys. 10, 864 (2014).
  * Steinhauer [2016] J. Steinhauer, Observation of quantum Hawking radiation and its entanglement in an analogue black hole, Nat. Phys. 12, 959–965 (2016).
  * Weinfurtner  _et al._ [2011] S. Weinfurtner, E. W. Tedford, M. C. J. Penrice, W. G. Unruh, and G. A. Lawrence, Measurement of stimulated Hawking emission in an analogue system, Phys. Rev. Lett. 106, 021302 (2011).
  * Visser and Weinfurtner [2005] M. Visser and S. Weinfurtner, Vortex analogue for the equatorial geometry of the Kerr black hole, Class. Quantum Gravity 22, 2493 (2005).
  * Torres  _et al._ [2017] T. Torres, S. Patrick, A. Coutant, M. Richartz, E. W. Tedford, and S. Weinfurtner, Rotational superradiant scattering in a vortex flow, Nature Physics 13, 10.1038/nphys4151 (2017).
  * Solnyshkov  _et al._ [2019] D. D. Solnyshkov, C. Leblanc, S. V. Koniakhin, O. Bleu, and G. Malpuech, Quantum analogue of a Kerr black hole and the Penrose effect in a Bose-Einstein condensate, Phys. Rev. B 99, 214511 (2019).
  * Braidotti  _et al._ [2022] M. C. Braidotti, R. Prizia, C. Maitland, F. Marino, A. Prain, I. Starshynov, N. Westerberg, E. M. Wright, and D. Faccio, Measurement of Penrose superradiance in a photon superfluid, Phys. Rev. Lett. 128, 013901 (2022).
  * Cromb  _et al._ [2020] M. Cromb, G. M. Gibson, E. Toninelli, M. J. Padgett, E. M. Wright, and D. Faccio, Amplification of waves from a rotating body, Nat. Phys. 16, 1069 (2020).
  * Jacquet  _et al._ [2023] M. J. Jacquet, L. Giacomelli, Q. Valnais, M. Joly, F. Claude, E. Giacobino, Q. Glorieux, I. Carusotto, and A. Bramati, Quantum vacuum excitation of a quasinormal mode in an analog model of black hole spacetime, Phys. Rev. Lett. 130, 111501 (2023).
  * Burgess  _et al._ [2024] C. Burgess, S. Patrick, T. Torres, R. Gregory, and F. König, Quasinormal modes of optical solitons, Phys. Rev. Lett. 132, 053802 (2024).
  * Patrick  _et al._ [2018] S. Patrick, A. Coutant, M. Richartz, and S. Weinfurtner, Black hole quasibound states from a draining bathtub vortex flow, Physical Review Letters 121, 10.1103/PhysRevLett.121.061101 (2018).
  * Torres  _et al._ [2020] T. Torres, S. Patrick, M. Richartz, and S. Weinfurtner, Quasinormal mode oscillations in an analogue black hole experiment, Physical Review Letters 125, 10.1103/PhysRevLett.125.011301 (2020).
  * Švančara  _et al._ [2024] P. Švančara, P. Smaniotto, L. Solidoro, J. F. MacDonald, S. Patrick, R. Gregory, C. F. Barenghi, and S. Weinfurtner, Rotating curved spacetime signatures from a giant quantum vortex, Nature 628, 66 (2024).
  * Torres  _et al._ [2022] T. Torres, S. Patrick, and R. Gregory, Imperfect draining vortex as analog extreme compact object, Phys. Rev. D 106, 045026 (2022).
  * Note [1] Note that our choice of $\varepsilon$ real and positive corresponds to a Neumann boundary condition in the case of a complete reflection, whereas choosing $\varepsilon$ negative/complex would lead to a Dirichlet/Robin boundary condition. In general, $\varepsilon$ could depend on $\omega$ but for the sake of simplicity, we model it here as a constant.
  * Macedo  _et al._ [2018] C. F. B. Macedo, T. Stratton, S. Dolan, and L. C. B. Crispino, Spectral lines of extreme compact objects, Phys. Rev. D 98, 104034 (2018).
  * [40] S. Dolan and T. Torres, Spectral (in)stability of black holes under potential perturbations, work in progress.
  * Florentin  _et al._ [1966] J. J. Florentin, M. Abramowitz, and I. A. Stegun, Handbook of mathematical functions., The American Mathematical Monthly 73, 10.2307/2314682 (1966).
  * Švančara  _et al._ [2024] P. Švančara, P. Smaniotto, L. Solidoro, J. F. MacDonald, S. Patrick, R. Gregory, C. F. Barenghi, and S. Weinfurtner, Rotating curved spacetime signatures from a giant quantum vortex, Nature , 1 (2024).
  * Barausse  _et al._ [2014] E. Barausse, V. Cardoso, and P. Pani, Can environmental effects spoil precision gravitational-wave astrophysics?, Phys. Rev. D 89, 104059 (2014).
  * Cole  _et al._ [2023] P. S. Cole, G. Bertone, A. Coogan, D. Gaggero, T. Karydas, B. J. Kavanagh, T. F. M. Spieksma, and G. M. Tomaselli, Distinguishing environmental effects on binary black hole gravitational waveforms, Nat. Astro. 7, 943 (2023).
  * Hannuksela  _et al._ [2019] O. A. Hannuksela, K. W. K. Wong, R. Brito, E. Berti, and T. G. F. Li, Probing the existence of ultralight bosons with a single gravitational-wave measurement, Nat. Astro. 3, 447 (2019).
  * Baumann  _et al._ [2022] D. Baumann, G. Bertone, J. Stout, and G. M. Tomaselli, Sharp signals of boson clouds in black hole binary inspirals, Phys. Rev. Lett. 128, 221102 (2022).
  * London  _et al._ [2014] L. London, D. Shoemaker, and J. Healy, Modeling ringdown: Beyond the fundamental quasinormal modes, Phys. Rev. D 90, 124032 (2014).
  * Mitman  _et al._ [2023] K. Mitman, M. Lagos, L. C. Stein, S. Ma, L. Hui, Y. Chen, N. Deppe, F. Hébert, L. E. Kidder, J. Moxon, _et al._ , Nonlinearities in black hole ringdowns, Phys. Rev. Lett. 130, 081402 (2023).

##  Appendix A Numerical simulations

We simulate (1) using the method of lines. First we define a spatial grid
$x_{k}$ with $k=1...N$ and $h=x_{k+1}-x_{k}$ and a 3-point centred difference
stencil for the second derivative,

$\partial_{x}^{2}\Psi_{k}=\frac{\Psi_{k+1}-2\Psi_{k}+\Psi_{k-1}}{h^{2}}.$ (12)

We then define $\Pi=\partial_{t}\Psi$ and cast (1) in the form,

$\partial_{t}\begin{pmatrix}\Psi\\\ \Pi\end{pmatrix}=\begin{pmatrix}0&1\\\
\partial_{x}^{2}-V&0\end{pmatrix}\begin{pmatrix}\Psi\\\ \Pi\end{pmatrix},$
(13)

and use a fourth-order Runge-Kutta algorithm to update the field at the next
time step. This method allows us to implement absorbing boundary conditions in
a simple way. At the edge of the domain where $V\simeq 0$ our solution will be
either a left-moving or right-moving plane wave, whose derivatives satisfy
$\partial_{x}\Psi=\pm\partial_{t}\Psi$. At the same order of accuracy in $h$,
the 3-point centred difference stencil for the first derivative,

$\partial_{x}\Psi_{k}=\frac{\Psi_{k+1}-\Psi_{k-1}}{2h},$ (14)

gives $\Psi_{0}=\Psi_{2}-2h\Pi_{1}$ at the left boundary and
$\Psi_{N+1}=\Psi_{N-1}-2h\Pi_{N}$ at the right boundary. Hence, we modify our
stencil for $\partial_{x}^{2}$ on the boundaries,

$\partial_{x}^{2}\Psi_{1}=\frac{2\Psi_{2}-2\Psi_{k-1}}{h^{2}}-\frac{2\Pi_{1}}{h}.$
(15)

and similarly for the point $k=N$. Note that without the term proportional to
$\Pi_{1}$, this would just be the second derivative implemented with Neumann
boundary conditions, say ${}^{\mathrm{N}}\partial_{x}^{2}$. Our numerical
equation of motion then gets modified to,

$\partial_{t}\begin{pmatrix}\Psi\\\
\Pi\end{pmatrix}=\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\
{}^{\mathrm{N}}\partial_{x}^{2}-V&B\end{pmatrix}\begin{pmatrix}\Psi\\\
\Pi\end{pmatrix},$ (16)

where $B=-(2/h)\mathrm{diag}(\mu,0,...0,1)$ is a diagonal matrix with non-zero
entries in the top-left and bottom-right corner. The parameter $\mu$ allows us
to tune between a fully reflective Neumann boundary condition on one side for
$\mu=0$ and a completely open boundary for $\mu=1$. The relation between $\mu$
and the reflectivity $\varepsilon$ used in the main text is,

$\mu=\frac{1-\varepsilon}{1+\varepsilon}.$ (17)

##  Appendix B Pöschl-Teller resonance condition

The linear wave equation (1) with the Pöschl-Teller potential reads:

$\displaystyle\frac{\partial^{2}\psi}{\partial
x^{2}}-\left(\omega^{2}-\frac{V_{0}}{\cosh(\alpha x)^{2}}\right)\psi=0\,.$
(18)

The equation can be reduced to an hypegeometric differential equation [41]:

$\displaystyle(1-y)y\frac{\partial^{2}\phi}{\partial
y^{2}}+[c-(a+b+1)y]\frac{\partial\phi}{\partial y}-ab\phi=0\,,$ (19)

whose solution is (Quasinormal modes in semi-open systems), by introducing the
variable:

$\phi=\psi\left((1-y)y\right)^{i\omega/2\alpha}\,.$

The parameters $a,b,c$ and the spatial variable $y$ are given in the main
text. In order to impose the radiating boundary condition (2) at
$x\to-\infty$, one sees that for large negative positions $y^{\pm
i\omega/2\alpha}\to e^{\pm i\omega x}$. Hence, it is possible to identify the
right-moving contribution in (Quasinormal modes in semi-open systems) and then
to impose the boundary condition on the left with $B=0$.

To look at the behaviour of the solution as $x\to+\infty$, i.e. $y\to 1$, one
can use the transformation rule of the hypergeometric functions [41] to get

$\psi=Ay^{-i\omega/2\alpha}\left[(1-y)^{i\omega/2\alpha}\frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}{}_{2}F_{1}(c-a,c-b,c^{*},1-y)+(1-y)^{-i\omega/2\alpha}\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}{}_{2}F_{1}(a,b,c,1-y)\right]\,,$

We can identify the out-going and in-going contributions as $x\to\infty$ by
noticing that

$(1-y)^{\pm i\omega/2\alpha}\xrightarrow{x\gg 1}e^{\mp i\omega x}\,.$

Therefore, one can define a Reflection Coefficient for the PT potential as the
ratio between the out-going and in-going contribution at the right of the
potential, resulting in (10)