删除或更新信息,请邮件至freekaoyan#163.com(#换成@)

The matrix method for black hole quasinormal modes

本站小编 Free考研考试/2022-01-01

Kai Lin 1,2,
, Wei-Liang Qian 2,3,4,
, 1.Hubei Subsurface Multi-scale Imaging Key Lab., Institute of Geophysics and Geomatics, China University of Geosciences, Hubei 430074, China
2.Escola de Engenharia de Lorena, Universidade de S?o Paulo, 12602-810, Lorena, SP, Brazil
3.Faculdade de Engenharia de Guaratinguetá, Universidade Estadual Paulista, 12516-410, Guaratinguetá, SP, Brazil
4.School of Physical Science and Technology, Yangzhou University, Jiangsu 225002, China
Received Date:2018-11-28
Available Online:2019-03-01
Abstract:We provide a comprehensive survey of possible applications of the matrix method for black hole quasinormal modes. The proposed algorithm can generally be applied to various background metrics, and in particular, it accommodates both analytic and numerical forms of the tortoise coordinates, as well as black hole spacetimes. We give a detailed account of different types of black hole metrics, master equations, and the corresponding boundary conditions. Besides, we argue that the method can readily be applied to cases where the master equation is a system of coupled equations. By adjusting the number of interpolation points, the present method provides a desirable degree of precision, in reasonable balance with its efficiency. The method is flexible and can easily be adopted to various distinct physical scenarios.

HTML

--> --> -->
1.Introduction
The black hole is one of the most exotic and intriguing subjects in all of theoretical physics. On the experimental side, black hole merger is a magnificent astrophysical phenomenon, it releases an enormous amount of energy in the form of gravitational radiation, larger even than that carried by the light from all stars in the entire visible Universe combined. The gravitational waves traverse the Universe, carrying with them the information about the merger. The first detection of gravitational waves in 2016 was heralded as a remarkable breakthrough in fundamental physics [1, 2]. The observation was announced by LIGO and Virgo collaborations, and it matches the predictions of gravitational waves that emanated from the inward spiral and merger of a pair of black holes [3-5]. Subsequently, further measurements confirmed the gravitational waves from compact-object binaries, including black holes and neutron stars [6-8]. These observations inaugurate a revolutionary era of astronomy, as astrophysicists and cosmologists are now able to make precise observations based on gravitational waves, in addition to electromagnetic radiation. Indeed, while black hole candidates have been identified through electromagnetic signals from nearby matter [9-12], gravitational wave measurements provide a direct test of general relativity, as well as unique access to the properties of spacetime, in the strong field regime. The final fate of a black hole binary consists of three distinct stages: inspiral, merger, and ringdown [13]. In the first stage of the inspiral, the orbit of the black hole binary shrinks gradually due to the emission of gravitational radiation. The inspiral dynamics can be calculated by using the post-Newtonian approximation, which results from a systematic expansion of the full Einstein equations. As the two black holes further evolve and spiral inward, they eventually reach the merger stage, which relates to the strong field, the dynamical regime of general relativity. During the merger stage, the two black holes coalesce into a single, highly distorted remnant black hole, which is surrounded by the combined event horizon. The process is no longer quasi-adiabatic and incredibly sophisticated, where most analytic or semi-analytic methods break down, and numerical relativity is usually employed to calculate the dynamical properties of the system. Gravitational wave emission reaches a peak during the second stage. The remnant black hole eventually settles down into a more dormant state, where distortions from the spherical shape rapidly decrease while emitting gravitational waves. This last stage is known as “ringdown”, in analogy with the way the ringing of a bell is suppressed in time after it is struck. The form of the gravitational wave is characterized by exponentially damped sinusoids. Various techniques in the black hole perturbation theory can be used to tackle the problem. The topic of the present paper are the black hole quasinormal modes, which are closely related to the last stage of the black hole merger process.
In general, a quasinormal mode is defined as an eigenmode of a dissipative system. In the context of general relativity, small perturbations of a black hole metric may lead to quasinormal modes [14-17]. Their importance is due to the fact that they are closely associated with the final fate of the black hole merger. Besides, quasinormal modes attracted much attention in recent years mainly due to the development of the holographic principle in the anti-de Sitter/conformal field theory (AdS/CFT) correspondence, which is widely recognized as a tool for exploring strongly coupled systems [18]. In other words, by studying the black hole perturbations in the bulk, one may extract important physical quantities, such as transport coefficients, e.g. viscosity, conductivity, and diffusion constants, of the dual system on the boundary.
The study of black hole quasinormal modes encompasses various types of perturbations, which include those of the scalar field, spinor field, vector field, or the metric itself. Revealedbytemporal evolution of small perturbations, the dynamic process of a stable black hole metric also involves three stages. These include the initial outburst, triggered by the source of the perturbation, the quasi-normal oscillations, and late-time asymptotic tail. Usually, the amplitude of the oscillations decays exponentially with time. Otherwise, if small oscillations increase, this indicates that the black hole metric is not stable. The frequency of a quasinormal mode is usually complex, and is known as quasinormal frequency. The real part of the quasinormal frequency defines the frequency of oscillations. The imaginary part, on the other hand, is negative for a stable metric, and it plays the role of suppressing the amplitude of oscillations.
From the mathematical point of view, the analysis of a quasinormal mode is related to a non-Hermitian eigenvalue problem for a system of coupled linear differential equations with associated boundary conditions. Apart of a few cases where analytic solutions for asymptotic quasinormal mode spectra can be obtained [19-21], one often resorts to semi-analytic techniques [22]. The P?shl-Theller potential approximation [23] replaces the effective potential in the master equation by a P?shl-Theller one, where the analytic eigenvalue is known. The WKB method [24, 25] is a semi-analytic approach proposed by Schutz et al, where the technique is applied specifically to the one-dimensional time-independent Schr?dinger-like equation for the quasi-eigenvalue problem. The 6th order formalism of the method was subsequently derived by Konoplya [26]. In practice, the application of the WKB method is rather straightforward, since it only requires the information about the background metric and the form of the effective potential. However, the uncertainty of the obtained result is somewhat undetermined beforehand, which is in part related to the specific shape of the effective potential. It is worth mentioning that, in the eikonal limit, the quasinormal modes of spherical black holes in asymptotically flat spacetime were found to be associated with the parameters of the circular null geodesic [27]. However, Konoplya and Stuchlik [28] showed analytically that the relationship is not valid for some particular cases, which has been further explored in [29-31]. In addition to the semi-analytic method, various numerical methods have also been proposed. The first type of numerical methods aims at evaluating the time-dependent wave equation directly for a given set of initial conditions. This is achieved by replacing the derivative with the finite difference [32-37]. Naturally, the resulting waveform contains the information about all three stages of temporal evolution of small perturbations. At the end of the calculation, one may further evaluate the quasinormal frequency by using Fourier analysis. A major disadvantage of the approach is that one cannot obtain the complete spectrum of quasinormal modes, due to the fact that usually only a few long-lived modes can be clearly identified. The second class of numerical schemes involves Fourier decomposition of the perturbed Einstein equations, while adopting appropriate boundary conditions. The latter assume that the resultant wave function corresponds to incoming waves on the horizon and outgoing at infinity. A series expansion of the wave function is usually carried out, and the method is typically iterative, namely, higher order results are derived from and improved using the lower order ones. Consequently, the precision of the method is controlled by the order of the expansion. These methods include the continued fraction method [38, 39], the Horowitz and Hubeny's (HH) method for the AdS black hole [40], and the asymptotic iteration method [41-44]. However, for the continued fraction method it is not so straightforward to derive the desired iterative relation for specific background metrics, and the HH method only applies to the asymptotically AdS spacetime. Recently, we proposed a matrix method [45-47] , where the spatial coordinate is discretized so that the differential equation, as well as its boundary conditions, are transformed into a homogeneous matrix equation. A vital feature of the method is that the eigenfunction is expanded in the vicinity of all grid points, and, therefore, the precision of the algorithm can be potentially improved. The method has been employed to study the quasinormal modes of the Schwarzschild black holes in asymptotically flat, (anti-)de Sitter (AdS/dS) spacetimes, as well as the Kerr and Kerr-Sen black hole metrics [46, 47]. In this survey, we further explore possible applications of the method to more general scenarios. We provide a detailed account of the procedure of how the algorithm is implemented. We argue that the implementation of the proposed method does not rely on the analytic form of the tortoise coordinate. Moreover, the method can be readily applied even to the case of the numerical black hole metric. Our discussion addresses black hole metrics of different asymptotic behavior, in particular where an analytic form of the tortoise coordinate is not available. Advantages, as well as improvements of the alternative approaches, are also addressed.
The paper is organized as follows. In the following section, we discuss the core algorithm of the matrix method. In Section III, we explore the applications of the method in various scenarios, including for different background metrics and the corresponding boundary conditions. In particular, we investigate the method of separation of variables for the master equation in asymptotically flat, dS, and AdS spacetimes, while considering different cases such as static, rotational, and extreme black holes. Moreover, the case where the master equation is a system of coupled ordinary differential equations is addressed. In addition, we provide a comprehensive discussion on the code implementation using Mathematica. Further discussions and concluding remarks are given in the last section of the paper.
2.The matrix method
The primary feature of the matrix method [45-47] , which makes it distinct from the others, lies in the fact that the series expansion of the wavefunction is carried out in the vicinity of a series of points, where the latter are not necessarily distributed evenly in the domain of the wavefunction. This feature provides flexibility for achieving a reasonable degree of precision, without sacrificing efficiency. Roughly speaking, in terms of the above grid point expansion, one discretizes the master equation of the perturbation field, and subsequently, rewrites it together with the boundary conditions in the form of a matrix equation. The latter can then be handled readily by various algorithms for solving a system of linear equations.
Let us start by discussing the general strategy of applying the method to the standard scenario for quasinormal modes, where the master equation is of the Schr?dinger-type, given as follows:
$\frac{{{{\rm d}^2}}}{{{\rm d}r_*^2}}\Phi (r) + \left[ {{\omega ^2} - V(r)} \right]\Phi (r) = 0,$
(2.1)
where $ r_*=\int{{\rm d}r/f(r)} $ is the tortoise coordinate, $ \omega $ and V(r) are the quasinormal frequency and effective potential, respectively. Here, f (r) is determined by the metric, while the black hole event horizon, rh , and the cosmological horizon, rc , are often determined by $ f(r_{\rm h})=0 $ and $ f(r_{\rm c})=0 $ , respectively. We proceed by noting that the above form of the master equation can generally be applied to a wide variety of metrics. However, it is not applicable neither to the case of rotating black holes, nor to that given in terms of Eddington–Finkelstein coordinates. The specific form of the master equation in these cases will be discussed with the corresponding topics in the following sections. Usually, the asymptotic behavior of the effective potential at the boundary is relatively simple, and, therefore, in the present work, we will just assume that the effective potential vanishes at the horizon and approaches a constant at infinity $ V(\infty)\equiv V_\infty=\rm{const} $.
In order to solve a differential equation, one also has to determine the associated boundary conditions. In the case of quasinormal modes, the boundary conditions are determined by the characteristics of the black hole at the horizon and infinity. As the apparent horizon of a black hole is a one-way membrane which prohibits the outgoing light rays to travel across it, the boundary conditions for the quasinormal modes are defined such that the solutions should be ingoing at the horizon and purely outgoing at infinity, namely
$\Phi \sim \left\{ {\begin{aligned} &{{{\rm e}^{{\rm i}\bar \omega {r_*}}}}&\!\!\!\!\!\!{r \to \infty }&\;\;{({\rm{for\ asymptotically \ flat \ spacetime}})}\\ &{{{\rm e}^{{\rm i}\omega {r_*}}}}&\!\!\!\!\!\!{r \to {r_{\rm c}}}&\;\;{({\rm{for \ asymptotically \ dS \ spacetime}})}\\ &{{{\rm e}^{ - {\rm i}\omega {r_*}}}}&\!\!\!\!\!\!{r \to {r_{\rm h}}}&\;\;{({\rm{for \ both \ spacetimes}})} \end{aligned}} \right.,$
(2.2)
where $ \bar{\omega} $ is defined as $ \bar{\omega}\equiv \sqrt{\omega^2-V_\infty} $, so that $ \bar{\omega}=\omega $ for $ V_\infty=0 $. If the form of the effective potential is not sufficiently simple, one usually has to resort to numerical methods. Eqs.(1) and (2) furnish an eigenvalue problem, where the boundary conditions are defined at $ r_*\rightarrow\pm\infty $ in terms of the tortoise coordinate. We note that the above discussion is valid for asymptotically flat and dS spacetimes, while the boundary conditions for asymptotically AdS spacetime, and the angular part of the master equation, will be addressed below.
Similarly to the other methods based on series expansion, we first change the domain of the radial variable to a finite range by introducing $ z=z(r) $. For convenience, we choose the new domain to be $ z\in[0,1] $, and the boundaries correspond to the points z=0 and z=1, respectively. However, since the wave function is also transformed, the boundary conditions in Eq.(2) will also be affected. In this regard, we rewrite the wave function in the form $ \Phi=A(z)R(z) $, where $ A(z(r)) $ carries the asymptotic behavior of the wave function at the boundary given in Eq.(2). As a result, R(z) is expected to be well behaved and, in particular, to be regular at the boundary. To be specific, the master equation reads
${\cal H}(\omega ,z)R(z) = 0,$
(2.3)
${\rm with}~ R(z = 0) = {C_0}~{\rm{and}}~ R(z = 1) = {C_1}.$
(2.4)
Here, $ {\cal H}(\omega,r) $ is a linear operator, determined by Eq.(1), which depends on $ \omega $ and z. C0 and C1 are constants, since the asymptotic part of the wave function has been factorized. It is convenient to further introduce
$F(z) = R(z)z(1 - z),$
(2.5)
where $ z(1-z) $ can be replaced by other smooth functions which are zero only at z=0 and z=1. We note that, for the case of AdS spacetime, since $ R(r\rightarrow\infty)\rightarrow 0 $, the transformation Eq.(5) can be replaced by $ F(z)=R(z)(1-z) $ where $ z=r_{\rm h}/r $.
Subsequently, the master equation can be rewritten with F(z)
${\cal G}(\omega ,z)F(z) = 0,$
(2.6)
with the boundary conditions in an even simpler form
$F(z = 0) = F(z = 1) = 0.$
(2.7)
Here, the linear operator $ {\cal G}(\omega,z) $ is derived from the form of $ {\cal H}(\omega,z) $ and Eq.(5). We note that, as a numerical trick, Eq.(5) is not an absolutely necessary procedure. In practice, the resultant wave function no longer possesses a boundary condition involving an undetermined constant. Furthermore, when C0 and C1 turn out to be quite distinct in their values, the factor $ {z(1-z)} $ helps to suppress the function values at the boundary points, as well as in their vicinity. Intuitively, the only scenario when the procedure does not work as intended is when the form of the wave function becomes more fluctuating due to the factor $ {z(1-z)} $for some particular reason, and the subsequent Taylor expansion becomes less efficient. But, as we do not have a good estimate of the form of the wave function at this point, the introduction of Eq.(5) helps to maintain the resultant equation simple, at the least. The existing numerical results are in most cases found to be consistent with the above speculations.
As the next step, we descretize the master equation Eq.(6) by introducing N interpolation points zi with $ i=1, 2, \cdots, N $ , in the interval $ z\in[0,1] $. Subsequently, the wave function F(z) is also represented by descrete values $ f_i=F(z_i) $ on the grid. As a result, one discretizes the differential equation Eq.(6) and rewrites it in terms of fi. In the present algorithm, various derivatives of the wave function at zi are obtained by inverting a $ N\times N $ matrix [45], so that the information about the function values in the entire interval is utilized. Besides, owing to the fact that the above inverting process is formal, it is carried out in practice before the real numerical calculations, and, therefore, it does not have any negative impact on the efficiency of the method. The resulting matrix equation can be formally written as
$ \bar{\cal M}{\cal F}=0, $
(2.8)
where $ \bar{\cal M} $ represents the descretized version of the operator $ {\cal G}(\omega,z) $, and is a $ N\times N $ matrix which is still a function of the quasinormal frequency $ \omega $, and
$ {\cal F}=(f_1,f_2,\cdots,f_i,\cdots,f_N)^T $
(2.9)
is a column vector composed of fi. The boundary condition Eq.(7) implies that
$ f_1=f_N=0. $
(2.10)
We make use of the above condition to replace the first and last line of the matrix $ \bar{\cal M} $. The reason behind this choice is to select a line or column in the original equation which makes use of the least amount of information about of grid points. Subsequently, Eq.(8) becomes
$ {\cal M}{\cal F}=0, $
(2.11)
where the element $ {\cal M}_{k,i} $ of the matrix $ {\cal M} $ is defined by
${{\cal M}_{k,i}} = \left\{ {\begin{array}{*{20}{l}} {{\delta _{k,i}},}&{k = 1\;{\rm{or}}\;{{N}}}\\ {{{\bar {\cal M} }_{k,i}},}&{k = 2,3, \cdots ,N - 1} \end{array}} \right..$
(2.12)
The matrix equation indicates that $ {\cal F} $ is the eigenvector of $ {\cal M} $, which implies
$ \det\left({\cal M}(\omega)\right)\equiv\left|{\cal M}((\omega)\right|=0. $
(2.13)
Eq.(13) is a non-linear algebraic equation for the quasinormal frequency $ \omega $, which can be solved by any standard nonlinear equation solver.
In Refs. [46, 47], the above algorithm is employed to descretize the master equation, Eq.(6). Alternatively, one may also employ other numerical approaches for numerical differentiation, such as the differential quadrature method [48-50], and the Runge-Kutta method, among others. Regarding the finite difference formulas, if Newton's difference quotient is employed in Eq.(6), the resulting matrix $ {\cal M}(\omega) $ will be “almost” diagonal, since only the elements $ {\cal M}_{k,i} $ with $ i=k, k\pm 1 $ are non-vanishing. Subsequently, the resulting wave function can be obtained via an iterative method, similar to the continued fraction method, or the HH method. However, as mentioned above, a sparse matrix indicates that the corresponding discretization scheme does not fully make use of the information about the entire interval $ z\in[0,1] $. On the other hand, in several different methods, such as the continued fraction method and the HH method, the series expansion of the wave function is only carried out in a single radial position, for instance, the horizon. From a different point of view, the resulting algebraic equation can in fact also be rephrased in the form of a matrix equation, Eq.(11). Owing to the fact that the iterative relation only involves a few expansion coefficients, the corresponding matrix is also rather sparse. Therefore, for a given N, the matrix method is related to a Taylor expansion of the order N, which subsequently provides better precision. Moreover, one may even seek out an optimized grid distribution, which provides a higher resolution for the region where the variation of the wavefunction is more dramatic, by inserting more grid points. In this regard, the present algorithm can be viewed as an improvement over many existing approaches.
The following section is devoted to a more detailed discussion of various applications.
3.Application of the method in the calculations of quasinormal frequency
While the general strategy for the evaluation of the black hole quasinormal modes presented in the previous section is concise and straightforward, many particular aspects, such as those related to the boundary conditions and specific characteristics of the master equations, might bring up complications for particular problems. In the present section, we discuss in detail several specific cases for the applications of the matrix method. In particular, our discussion addresses the boundary conditions for the black hole metrics in asymptotically flat, dS, and AdS spacetimes, as well as the associated ansatz for the method of separation of variables. The background metrics concern static, rotational, and extreme black holes. The algorithm presented in the following does not require the existence of an analytic form of the tortoise coordinate. We also comment on the case of a numerical black hole metric, and that where the master equation is a system of coupled ordinary differential equations.
We first address the issue of different boundary conditions owing to different asymptotical properties of spacetime. Usually, spacetimes can be classified according to their asymptotical behavior at infinity, such as asymptotically flat, dS, and AdS spacetimes. It is known that the effective potential V(r) vanishes at infinity in the case of asymptotically flat spacetime. On the other hand, it diverges at infinity for asymptotically dS and AdS spacetimes. However, for the dS case, there exists a cosmological horizon rc between the observer and infinity, which is also a Cauchy one-way membrane. An observer staying between rh and rc will not receive any signals from the inside of event horizon, nor from the outside of cosmological horizon. Therefore, the physically valid domain is $ r_{\rm h}\leqslant r \leqslant r_{\rm c} $ for the dS spacetime, while it is $ r_{\rm h} \leqslant r \leqslant \infty $ for the asymptotically flat and AdS spacetimes.
For extreme black holes, the associated tortoise coordinate has a different feature near the black hole horizon than for non-extreme black holes. This is because, for an extreme black hole, one requires $ f'(r_{\rm h})=0 $, which implies a Taylor expansion of f (r) around the horizon that does not contain the linear term. Subsequently, as shown below, the leading contribution in the tortoise coordinate becomes different.
For rotating black holes, an additional equation is implied, whose physical content is associated with the angular momentum quantum number. In this case, one usually gets a system of two coupled master equations with eigenvalues $ \omega $ and $ \lambda $, corresponding to the radial and angular parts of the wave functions. The latter restores to a simple case when the parameter of the black hole related to angular asymmetry vanishes. In particular, if the black hole metric reduces to a static spherical one, one finds $ \lambda=\ell(\ell+1) $, namely, the angular momentum quantum number. The angular part of the master equation becomes decoupled, and the angular wavefunction reduces to spherical harmonics. On the other hand, if the angular part of the master equation is coupled to the radial part, its boundary condition is the usual periodic boundary condition.
In addition, one frequently encounters two difficult scenarios in realistic numerical studies. The first involves the difficulty in obtaining an analytic expression for the tortoise coordinate $ r_*(r) $ , and as a result it is not straightforward to obtain the boundary condition in an analytic form. The second difficulty is apparently even more severe, and it concerns the numerical black hole metric. Fortunately, as we show below, the matrix method is quite versatile in handling these subtle situations. In fact, it is relatively intuitive to understand the reason. The critical procedure of the matrix method is the discretization of continuous functions into discrete values on a grid. Naturally, as a result, the implementation of the method does not rely on analytic expressions of the tortoise coordinate, as shown in the following subsection. Moreover, the method neither relies on the analytic form of the black hole metric itself. It is rather straightforward to slightly alter the procedure to adopt it to the numerical black hole metric. This is not the case with the other approaches, such as the HH method, where the numerical difficulty is enhanced by the requirement of precision values for the high order derivatives of the metric. Last but not least, for the difficulties with the numerical tortoise coordinate, one does not require an explicit form of r* but its asymptotical behavior, and, therefore, one can merely expand the form of r* near the boundary and extract the desired boundary conditions.
2
3.1.Boundary conditions
-->

3.1.Boundary conditions

In the following, we write down the explicit forms of the boundary conditions for five specific cases of asymptotic spacetime and its boundaries:
Case 1. For the non-extreme black hole, near the horizon, if $ f(r)=\hat{f}_1(r-r_{\rm h})+\displaystyle\frac{\hat{f}_2}{2}(r-r_{\rm h})^2+{\cal O}(r-r_{\rm h}) $ (where $ \hat{f}_i\equiv f^{(i)}(r_{\rm h}) $ and $ {\cal O} $ represents terms of order higher than $ (r-r_{\rm h})^2 $), we have
$\begin{split} {r_*} = &\int {\displaystyle\frac{{{\rm d}r}}{{f(r)}}} = \int {\left( {\displaystyle\frac{1}{{{{\hat f}_1}(r - {r_{\rm h}})}} - \frac{{{{\hat f}_2}}}{{2\hat f_1^2}} + O} \right)} {\rm d}r\\ =& \displaystyle\frac{{\ln \left| {r - {r_{\rm h}}} \right|}}{{{{\hat f}_1}}} - \displaystyle\frac{{{{\hat f}_2}}}{{2\hat f_1^2}}(r - {r_{\rm h}}) + O, \end{split}$
(3.1)
which implies that the boundary condition at the horizon is given by
$ \Phi\sim {\rm e}^{-{\rm i}\omega r_*}\sim B_a\equiv \left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}}. $
(3.2)
Case 2. For the extreme black hole, near the horizon, if $ f(r)\!=\!\displaystyle\frac{\hat{f}_2}{2}(r\!-\!r_{\rm h})^2\!+\!\displaystyle\frac{\hat{f}_3}{6}(r\!-\!r_{\rm h})^3\!+\!\displaystyle\frac{\hat{f}_4}{24}(r\!-\!r_{\rm h})^4\!+\!{\cal O}(r\!-\!r_{\rm h}) $, where $ {\cal O} $ are terms of order higher than $ (r-r_{\rm h})^4 $, we find
$\begin{split} {r_*} &= \int {\displaystyle\frac{{{\rm d}r}}{{f(r)}}} \\&= \int {\left( {\displaystyle\frac{2}{{{{\hat f}_2}{{(r - {r_{\rm h}})}^2}}} - \displaystyle\frac{{2{{\hat f}_3}}}{{3\hat f_2^2(r - {r_{\rm h}})}} + \displaystyle\frac{{4\hat f_3^2 - 3{{\hat f}_2}{{\hat f}_4}}}{{18\hat f_2^3}} + O} \right)} {\rm d}r\\ &= - \displaystyle\frac{2}{{{{\hat f}_2}(r - {r_{\rm h}})}} - \displaystyle\frac{{2{{\hat f}_3}\ln \left| {r - {r_{\rm h}}} \right|}}{{3\hat f_2^2}} + \displaystyle\frac{{4\hat f_3^2 - 3{{\hat f}_2}{{\hat f}_4}}}{{18\hat f_2^3}}(r - {r_{\rm h}}) + O, \end{split}$
(3.3)
and the boundary condition at horizon is given by
$ \Phi\sim {\rm e}^{-{\rm i}\omega r_*}\sim B_b\equiv\left(r-r_{\rm h}\right)^{\frac{2i\omega \hat{f}_3}{3\hat{f}_2^2}}{\rm e}^{\frac{2i\omega}{\hat{f}_2(r-r_{\rm h})}}. $
(3.4)
Case 3. In the asymptotically flat spacetime, for $ r\to \infty $, if $ f(r)=1+A r^\alpha+{\cal O} $, where $ \alpha<0 $ and $ {\cal O} $ represents higher order terms less significant than $ r^{\alpha} $, we have
${r_*} = \int {\frac{{{\rm d}r}}{{f(r)}}} = {\left. {\int {\left( {1 - A{r^\alpha } + O} \right)} {\rm d}r} \right|^{r \to \infty }} = r - \frac{A}{{1 + \alpha }}{r^{1 + \alpha }} + O.$
(3.5)
Thus, the boundary condition at infinity is given by
$ \Phi\sim {\rm e}^{{\rm i}\bar{\omega} r_*}\sim B_0\equiv {\rm e}^{{\rm i}\bar{\omega}\left(r-\frac{A}{1+\alpha}r^{1+\alpha}\right)},$
(3.6)
and the second term could be significant when $ \alpha>-1 $.
Case 4. In the asymptotically de Sitter spacetime, according to Case 1 and Case 2, one obtains similar results by replacing rh with rc. Therefore, following Eq.(15), near the non-extreme cosmological horizon, the boundary condition is
$ \Phi\sim {\rm e}^{i\omega r_*}\sim B_1\equiv\left(r_{\rm c}-r\right)^{\frac{i\omega}{\bar{f}_1}}. $
(3.7)
On the other hand, for the extreme cosmological horizon, we have
$ \Phi\sim {\rm e}^{{\rm i}\omega r_*}\sim B_2\equiv\left(r_{\rm c}-r\right)^{-\frac{2i\omega \bar{f}_3}{3\bar{f}_2^2}}{\rm e}^{\frac{2{\rm i}\omega}{\bar{f}_2(r_{\rm c}-r)}}, $
(3.8)
where $ \bar{f}_i\equiv f^{(i)}(r_{\rm c}) $.
2
3.2.Separation of variables
-->

3.2.Separation of variables

By making use of the asymptotic behavior obtained above, we can write down the appropriate forms for the wave functions by the method of separation of variables.
For the non-extreme black hole in asymptotically flat spacetime, we have
$ \Phi(r)={\rm e}^{{\rm i}\bar{\omega}\left(r-\frac{A}{1+\alpha}r^{1+\alpha}\right)}\left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}}r^{\frac{i\omega}{\hat{f}_1}}R(r), $
(3.9)
where the term $ r^{\frac{i\omega}{\hat{f}_1}} $ is introduced to cancel out the effect of the factor $ \left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}} $ at infinity, while it gives a rather smooth contribution at the horizon. Thus, it can be readily verified that the resultant expression, Eq.(22), indeed possesses the desired asymptotic forms at $ r\to r_{\rm h} $ and $ r\to\infty $ found in Eqs.(15) and (18). Likewise, for the extreme black hole, we have
$\Phi (r) = {{\rm e}^{{\rm i}\bar \omega \left( {r - \frac{A}{{1 + \alpha }}{r^{1 + \alpha }}} \right)}}{\left( {r - {r_{\rm h}}} \right)^{\frac{{2i\omega {{\hat f}_3}}}{{3\hat f_2^2}}}}{{\rm e}^{\frac{{2{\rm i}\omega }}{{{{\hat f}_2}(r - {r_{\rm h}})}}}}{r^{ - \frac{{2i\omega {{\hat f}_3}}}{{3\hat f_2^2}}}}{{\rm e}^{ - \frac{{2{\rm i}\omega }}{{{{\hat f}_2}r}}}}R(r).$
(3.10)
As discussed above, for the black hole in asymptotically de Sitter spacetime, it is required that at the event horizon $ \hat{f}_1\not=0 $ and $ \hat{f}_1=0 $ for the non-extreme and extreme cases, respectively. Similarly, at the cosmological horizon, we have $ \bar{f}_1\not=0 $ and $ \bar{f}_1=0 $ , respectively for the non-extreme and extreme scenarios. Therefore, the boundary conditions in de Sitter spacetime are satisfied by assuming the following form for the wave functions
$\Phi (r) = \left\{ {\begin{array}{*{20}{c}} {{B_a}{B_1}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{\not = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{\not = 0}}}\\ {{B_a}{B_2}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{\not = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{ = 0}}}\\ {{B_b}{B_1}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{ = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{\not = 0}}}\\ {{B_b}{B_2}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{ = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{ = 0}}} \end{array}} \right..$
(3.11)
For the asymptotically anti-de Sitter spacetime, the Horowitz-Hubeny method usually employs the ingoing Eddington-Finkelstein coordinates $ (v,r) $ , with $ v=r_*+t $. In what follows, we will also use this convention. The corresponding master equation reads
$ f(r)\frac{{\rm d}^2\Phi(r)}{{\rm d}r^2}+\left[\frac{{\rm d}f(r)}{{\rm d}r}-2i\omega\right]\frac{{\rm d}\Phi(r)}{{\rm d}r}-U(r)\Phi(r)=0, $
(3.12)
where $ U(r)=V(r)/f(r) $, and since the potential is divergent at infinity, the boundary condition requires $ \Phi $ to be constant at the event horizon, while it must vanish at infinity.
For the rotating black hole, one has to deal with the angular part of the master equation, which possesses the following form, similar to the radial part of the master equation:
$ (1-u^2)\frac{{\rm d}}{{\rm d}u}\left[(1-u^2)\frac{{\rm d}Y(u)}{{\rm d}u}\right]-{\cal W}(u,\omega,\lambda)Y(u)=0. $
(3.13)
The boundary conditions are at $ u=\pm1 $ , where one defines $ W^2_\pm ={\cal W}(u=\pm 1,\omega,\lambda) $. The asymptotic solutions at the boundary are given by
$Y(u) \sim \left\{ {\begin{array}{*{20}{l}} {{{\left| {1 - u} \right|}^{ \pm \frac{{{W_ + }}}{2}}}}&{u \to 1}\\ {{{\left| {1 + u} \right|}^{ \pm \frac{{{W_ - }}}{2}}}}&{u \to - 1} \end{array}} \right..$
(3.14)
As it is required that the angular part of the wave function Y(u) always remains finite, it can be assumed to have the following form
$Y(u) = {\left| {1 - u} \right|^{\left| {\frac{{{W_ + }}}{2}} \right|}}{\left| {1 + u} \right|^{\left| {\frac{{{W_ - }}}{2}} \right|}}S(u).$
(3.15)
In the above expressions for the wave function with separation of variables, Eqs.(22)-(28), we have implied that the functions R(r) and S(u) are finite and well behaved in the entire domain of r and u. That is to say, $ r_{\rm h}\leqslant r < \infty $ for the asymptotically flat and anti-de Sitter spacetimes, $ r_{\rm h}\leqslant r\leqslant r_{\rm c} $ for the asymptotically de Sitter spacetime, and $ -1\leqslant u\leqslant 1 $ for the angular part of the master equation. This is because the asymptotic form of the wave function at the boundary has already been factorized and, therefore, effectively extracted from the l.h.s. of Eqs.(22)-(28).
Last but not least, we carry out another coordinate transformation to conveniently convert the domain of r and u into [0,1]. For the asymptotically flat and anti-de Sitter black hole spacetimes, the new coordinate is chosen as $ x=1-\displaystyle\frac{r_{\rm h}}{r} $ or $ z=r_{\rm h}/r $, where rh is the event horizon, so that at the event horizon x=0 and z=1, while x=1 and z=0 represent infinity. For the asymptotically de Sitter black hole spacetime, due to the existence of the cosmological horizon rc, we choose $ y=\displaystyle\frac{r-r_{\rm h}}{r_{\rm c}-r_{\rm h}} $ , so that y=0 represents the event horizon while y=1 corresponds to the cosmological horizon. For the angular coordinate, we introduce $ \nu=\displaystyle\frac{1+u}{2} $ , so that one has $ \nu=0 $ and $ \nu=1 $ for u=?1 and u=1, respectively.
By following the above procedure, the transformed master equation and its boundary conditions are given in Eq.(3) and Eq.(4). Subsequently, one may proceed with the calculations of the quasinormal frequency of the matrix equation as discussed in Section 2.
2
3.3.Master equation for a system of coupled equations
-->

3.3.Master equation for a system of coupled equations

In the above discussion, we have assumed that the master equation can be expressed in terms of a system of decoupled equations, as we have treated the radial and angular parts of the master equation separately. In practice, such simplification might not always be possible. As an example, in the study of the quasinormal modes of a massive vector field, the resultant system of equations might be coupled. However, a system of coupled ordinary differential equations does not particularly pose a difficulty for the matrix method. This is because, in this case, one only needs to write the system of coupled equations in terms of a matrix equation of larger size. In other words, if a system of n equations cannot be decoupled, the complication formally manifests in the fact that one has to handle a matrix equation of the size $ nN\times nN $, where N is the number of interpolation points introduced above. As an example, in the case of n=2, we have
$\begin{array}{l} {{\cal M}_a}{{\cal F}_a} = {{\cal B}_b}{{\cal F}_b},\quad {{\cal M}_b}{{\cal F}_b} = {{\cal B}_a}{{\cal F}_a}, \end{array}$
(3.16)
where $ {\cal B} $ is also a $ N\times N $ matrix. In this case, one may rewrite the above equations in terms of a double-sized matrix equation, namely,
${{\cal M}_n}{{\cal F}_n} \equiv \left( {\begin{array}{*{20}{c}} {{{\cal M}_a}}&{ - {{\cal B}_b}}\\ { - {{\cal B}_a}}&{{{\cal M}_b}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\cal F}_a}}\\ {{{\cal F}_b}} \end{array}} \right) = 0,$
(3.17)
which can be readily solved as before. The only disadvantage is that the matrix on the l.h.s. of Eq.(30) might be a sparse matrix, which subsequently implies a loss of computational efficiency.
Alternatively, one can either rewrite Eq.(30) as
$ \left({\cal M}_a{\cal B}_a^{-1}{\cal M}_b-{\cal B}_b\right){\cal F}_b=0, $
(3.18)
or as
$ \left({\cal M}_b{\cal B}_b^{-1}{\cal M}_a-{\cal B}_a\right){\cal F}_a=0. $
(3.19)
Here, the expressions involve the evaluation of $ {\cal B}_a^{-1} $ or $ {\cal B}_b^{-1} $, when the relevant inverse exists. To solve Eq.(31) or Eq.(32) , it is necessary to find the roots of a nonlinear equation in $ \omega $, associated with the $ N\times N $ matrices. The second approach may potentially improve the efficiency.
2
3.4.Code implementation
-->

3.4.Code implementation

The code implementation has been carried out in Mathematica. It takes advantage of the efficiency of existing packages for matrix manipulation, as well as for solving nonlinear equations. The source codes for the Schwarzschild black holes in asymptotically flat, (anti-) de Sitter (AdS/dS) spacetimes, as well as for the Kerr and Kerr-Sen black hole metrics, have been released publicly [46, 47]. In this subsection, we discuss few technical aspects of the codes implemented in Mathematica.
Firstly, Mathematica is not always very efficient in handling decimals, especially when decimals are further processed as arguments of internal functions, such as trigonometric or exponential functions. Therefore, in view of optimizing the computational time, one should avoid coding in decimals. As an example, one may use the following Macro to transform decimal output of a function “TU” into a fraction:
$ \begin{split} \rm{YOULI}[\rm{TU}\_] :=& \rm{Rationalize}[\rm{Chop}[\rm{Expand}\\ &[\rm{N}[\rm{TU},\rm{precs}]],10^{-\rm{precs}}], 10^{-\rm{precs}}], \end{split} $
(3.20)
where “precs” is the desired precision of the function “TU”. If one assigns precs=100, the function output will be given in terms of a fraction truncated to a hundred significant digits.
In Mathematica, there are usually two methods of finding the roots of Eq.(13), which is essentially an algebraic equation. They are FindRoot and NSolve. The former is used to search for a root of a function starting from an initial point provided by the user, while the latter attempts to enumerate all existing roots of the function. For the first case, one has to find an appropriate initial guess for the quasinormal frequency, which can be achieved by the following approaches. (1) One can make use of the result of another less accurate method, such as P?shl-Theller potential approximation, as input for the FindRoot command. (2) When the black hole metric in question restores to a more straightforward case, for which the quasinormal frequency is already known, by continuously varying its parameters, one may obtain the desired result by gradually changing the parameters of the metric to the relevant values. We note that the above approaches may not work as intended for the case where the quasinormal frequency possesses a bifurcation in the parameter space. However, this problem exists in general for a variety of methods where the eigenvalue problem is expressed in terms of nonlinear equations. (3) One may also make use of the NSolve package to calculate the quasinormal frequencies. In principle, one may find quasinormal frequencies related to different quantum numbers by a single root-finding process. However, the challenge is how to pick out real physical roots in the numerical results. One should proceed with care, and compare the numerical values with those obtained from known cases of black hole metrics. Moreover, the following rule of thumb may provide some general guidance. For some cases, in particular those where the resultant algebraic equation, Eq.(13), is merely an N-th order polynomial, one can exhaust all roots via the NSolve command. Among these roots, some are not physical, and they appear because of a certain symmetry of the equations associated with a particular choice of interpolation points. These unphysical roots can be detected by varying the number of interpolation points, since the roots corresponding to the physical quasinormal frequency will always be present and approach a limit as $ N\to \infty $. Nonetheless, to properly remove all unphysical roots is not a simple task. The situation becomes even more complicated when Eq.(13) is highly nonlinear. In this regard, the above mentioned Macro can be employed to alleviate the situation.
Apart from the above two methods, the matrix method can also be adapted to an iterative scheme. This third method is meaningful especially for situations where neither FindRoot nor NSolve is found to be efficient. In general, the iterative relation can be written as follows.
$ \omega_{k+1}={\cal Y}\left(\omega_k,\left|{\cal M}(\omega)\right|\right), $
(3.21)
where the quasinormal frequency of the (k+1)-th interation, $ \omega_{k+1} $, is determined by the k-th result. The operator $ {\cal Y} $ is a functional of the determinant satisfying $ \omega={\cal Y}\left(\omega,\left|{\cal M}\right|\right) $ , where $ \omega $ is the root of Eq.(13). By choosing an appropriate form of $ {\cal Y} $, one obtains the result for quasinormal frequency by repeatedly using Eq.(25). The primary challenge of this approach is to find a proper iteration relation for $ {\cal Y}(x,y) $ which efficiently leads to a convergent result.
In our code implementation in Refs.[46, 47], we have evenly divided the domain [0,1] ofz. In practice, it might be more favorable to utilize a nonuniform distribution which introduces more interpolation points in the region where the effective potential changes more radically. As a matter of fact, this is an issue that one encounters in other numerical schemes, such as the differential quadrature method and the Runge-Kutta method.
2
3.5.Additional considerations
-->

3.5.Additional considerations

There are a few more issues or topics which have not been explored in the existing studies but are worth mentioning.
Apart for the asymptotically AdS spacetime, our analysis has been mostly carried out in the radial coordinate r. However, sometimes it can be more convenient to express the formalism in the ingoing Eddington–Finkelstein coordinates. In this case, the corresponding boundary condition has also to be modified
$ \Phi\sim \left\{\begin{array}{*{20}{l}} {\rm e}^{{\rm i}\left(\omega+\tilde{\omega}\right) r_*} & r\rightarrow\infty ~\,{\rm{or}}\,~ r\rightarrow r_{\rm c}\\ \Phi_0 & r\rightarrow r_{\rm h} \end{array}\right.. $
(3.22)
where $ \tilde{\omega}=\sqrt{\omega^2-V_\infty} $ for $ r\rightarrow\infty $ in the asymptotically flat spacetime, and $ \tilde{\omega}=\omega $ for $ r\rightarrow r_{\rm c} $ in the asymptotically dS spacetime. We note that, in a similar fashion, the master equation can also be derived in the outgoing Eddington-Finkelstein coordinates.
In practice, instead of carrying out the transformation Eq.(5) and rewriting the boundary condition as Eq.(7), the original boundary condition Eq.(4) might be directly employed in the calculations. As mentioned before, the primary motivation of introducing Eq.(5) is when C0 and C1 are very distinct in their magnitudes. It is understood that such a procedure is useful for the radial coordinate r , as discussed in the present work.
In reality, the black hole is not a static object. The topic of quasinormal modes for dynamic black holes is not only intriguing but also significant in practice. For a dynamical black hole, the event horizon, apparent horizon and infinite redshift surface are distinct concepts, and, therefore, the boundary conditions for quasinormal modes should be tackled with special care.
In certain modified gravity theories, the Lorentz invariance is broken. Owing to the tachyon, the concepts of the event horizon and apparent horizon have to be modified in these theories. Recent studies have proven the existence of the universal horizon [51, 52], which is defined as the boundary that even superluminal particles cannot escape to future infinity. As a part of the effort to find evidence of Lorentz invariance breaking, it is, therefore, interesting to investigate the quasinormal modes for the universal horizon in the framework of the matrix method.
The notion of quasinormal modes is not restricted to black holes. Quasinormal modes of compact stars have since long aroused considerable attention. The primary difference between the two phenomena, from the mathematical point of view, is the existence of the outgoing wave at the boundary of the star surface. This implies a more subtle condition for properly connecting the inside and outside wave functions at the boundary. Nonetheless, the matrix method can also be adapted to handle such problems appropriately.
Last but not least, we note that a quasinormal mode by itself is an eigenvalue of a secular equation. The matrix method is merely a discretized version of the eigenvalue problem in question. One of the alternatives to the eigenvalue problem is the shooting method, which solves the boundary value problem by reducing it to the initial value problem. In this context, one might also implement a specific version of the shooting method for the column vector, which is nothing else but a discretized eigenfunction.
4.Concluding remarks
In this paper, a comprehensive survey of the applications of the matrix method for black hole quasinormal modes was presented. It was shown that the proposed algorithm is independent of any particular analytical form of the tortoise coordinate. In fact, the method can even be employed to a numerical black hole metric. Moreover, it is argued that the present approach can readily be applied to the case where the master equation involves a system of coupled equations. Our discussion covered various types of black hole metrics, master equations, and the corresponding boundary conditions. The proposed method is reasonably accurate, efficient, as well as flexible for distinct physical scenarios. We look forward to carrying out further studies of several open questions addressed in this survey.
闂傚倸鍊搁崐鐑芥嚄閸撲礁鍨濇い鏍仜缁€澶嬩繆閵堝懏鍣圭紒鐘靛█閺岀喖骞戦幇闈涙闂佸憡淇洪~澶愬Φ閸曨垰妫橀柛顭戝枓閹稿啴姊洪崨濠庢畷鐎光偓閹间礁绠栨俊銈呮噺閺呮煡骞栫€涙ḿ绠橀柣鈺佹捣缁辨挻鎷呮搴ょ獥闂侀潻缍囩紞浣割嚕婵犳碍鍋勯柣鎾虫捣椤ρ囨⒑閸忚偐銈撮柡鍛箞閹繝宕掗悙绮规嫼缂備礁顑堝▔鏇㈡倿閸ф鐓欓柛鎴欏€栫€氾拷2婵犵數濮烽弫鎼佸磻閻愬搫鍨傞柛顐f礀缁犳壆绱掔€n偓绱╂繛宸簻鎯熼梺鍐叉惈椤戝洨绮欒箛娑欌拺闁革富鍘奸崝瀣亜閵娿儲顥㈢€规洜鏁婚崺鈧い鎺戝閳锋垿鏌涘☉姗堝伐濠殿噯绠戦湁婵犲﹤鎳庢禒杈┾偓瑙勬礃濡炰粙寮幘缁樺亹鐎规洖娲ら獮妤呮⒒娓氣偓濞佳呮崲閸儱纾归柡宓偓濡插牏鎲搁弮鍫濊摕闁挎繂顦悞娲煕閹板吀绨奸柛锝勫嵆濮婅櫣鎷犻垾铏闂佹悶鍎滈崶褎鏆梻鍌欑劍鐎笛呮崲閸屾娲閵堝懐锛涢梺鍦劋椤ㄥ棝鍩涢幋锔界厱婵犻潧妫楅鈺呮煃瑜滈崜娆戠礊婵犲洤绠栭梺鍨儐缂嶅洭鏌嶉崫鍕簽婵炶偐鍠栧铏规崉閵娿儲鐝㈤梺鐟板殩閹凤拷
婵犵數濮烽弫鍛婃叏娴兼潙鍨傜憸鐗堝笚閸嬪鏌曡箛瀣偓鏇㈡倷婵犲嫭鍠愮€广儱妫欓崣蹇涙煏閸繍妲归柍閿嬪灴閺屾稑鈽夊鍫濅紣缂備焦顨嗙敮妤佺┍婵犲浂鏁冮柨婵嗘处閸掓稑顪冮妶鍐ㄧ仾婵☆偄鍟幈銊╁焵椤掑嫭鐓忛柛顐g箖閿涘秵淇婇銏狀伃闁哄矉绲鹃幆鏃堫敍濠婂憛锝夋⒑閸濄儱校闁绘濮撮悾鐑藉閵堝懐顔掑銈嗘⒒閺咁偊宕㈤幖浣光拺闁告稑锕ョ粈瀣箾娴e啿娲﹂崐鍫曟煥濠靛棙顥撳ù婊勭矒閺岀喓鈧稒岣跨粻鏍ь熆鐠哄搫顏紒杈ㄥ笧閳ь剨缍嗘禍璺何熼埀顒勬⒑缁洘鏉归柛瀣尭椤啴濡堕崱妤€娼戦梺绋款儐閹瑰洭寮诲鍥ㄥ珰闁哄被鍎卞鏉库攽閿熺姷鐣哄ù婊冪埣瀵顓奸崼顐n€囬梻浣告啞閹稿鎮烽埡浣烘殾妞ゆ牗绋戦閬嶆倵濞戞顏呯椤栨埃鏀介柣鎰级閳绘洖霉濠婂嫮绠炵€殿喗鐓¢、妤呭礋椤掆偓閳ь剙鐖奸弻锝夊箛椤旇姤姣勯梺纭呮閸婂潡寮诲☉銏犖ч柛銉仢閵忋倖顥嗗璺侯儑缁♀偓婵犵數濮撮崐鎼佸汲閿濆棎浜滈幖娣焺濞堟洟鏌曢崶褍顏柛鈺冨仱椤㈡﹢鎮欏顔荤棯濠电姵顔栭崹閬嶅箰閹惰棄钃熼柨鐔哄Т閻愬﹪鏌嶆潪鎵妽闁诲繋绶氬娲川婵犲嫭鍠涢梺绋款儐閹瑰洤顫忕紒妯诲闁告縿鍎虫婵犵數鍋橀崠鐘诲幢閹邦亝鐫忛梻浣虹帛閸旀寮崫銉т笉闁哄啫鐗婇悡娆撴煙椤栧棗鑻▓鍫曟⒑瀹曞洨甯涙慨濠傜秺楠炲牓濡搁妷搴e枔閹风娀骞撻幒婵囨祰闂傚倷鐒﹂幃鍫曞磹瑜忕划濠氬箻鐠囪尪鎽曢梺缁樻濞咃綁鎯屽▎鎾寸厵缂佸鐏濋銏ゆ煙椤旂晫鎳囨慨濠勫劋鐎电厧鈻庨幋鐘樻粎绱撴担鍝勑i柣妤佹礋椤㈡岸鏁愭径妯绘櫇闂佸啿鐏堥弲婊堟倵婵犳碍鈷戠憸鐗堝笒娴滀即鏌涘Ο鍝勨挃缂侇喗鐟╁畷鐔碱敍濞戞帗瀚奸梻浣告贡鏋繛瀵稿厴閸┿儲寰勯幇顓犲幐闂佸壊鍋掗崑鍕櫠鐎电硶鍋撶憴鍕缂傚秴锕ユ穱濠傤潰瀹€濠冃┑鐘愁問閸ㄤ即濡堕幖浣歌摕婵炴垶菤濡插牊鎱ㄥΔ鈧悧濠囧极閸撗呯=濞达絽鎼牎闁汇埄鍨抽崑銈夊春閳ь剚銇勯幒鍡椾壕闂佽绻戦懝楣冣€﹂崹顕呮建闁逞屽墴楠炲啳顦圭€规洖宕湁闁哄瀵ч崰妯尖偓瑙勬礈鏋摶鏍归敐澶嬫珳闁汇儺浜缁樻媴娓氼垱鏁梺瑙勬た娴滎亜顫忔禒瀣妞ゆ牗绋掑▍鏍⒑閸濆嫮鈻夐柛妯圭矙閹ょ疀濞戞瑧鍘遍梺鏂ユ櫅閸燁垳绮堥埀顒€顪冮妶蹇曞矝闁哄棙绔糴婵犵數濮烽弫鍛婃叏娴兼潙鍨傞柛锔诲幘缁€濠傗攽閻樺弶鎼愰柣鎺戠仛閵囧嫰骞掑鍫濆帯闂佹剚鍨卞ú鐔煎蓟閺囥垹骞㈡俊銈傚亾闁哄棴缍侀弻锛勪沪閸撗勫垱濡ょ姷鍋炵敮锟犵嵁鐎n喗鍊婚柛鈩冿供濡冣攽閿涘嫬浜奸柛濠冪墱閺侇噣鎮欓崫鍕崶闂佸綊鍋婇崰姘舵儗濞嗗繆鏀介柣妯哄级婢跺嫰鏌涚€n偄濮嶉柡宀嬬秮婵偓闁靛繆鍓濆В鍕煛娴e摜澧︽慨濠勭帛閹峰懐绮欓幐搴♀偓顖氣攽閻橆喖鐏柨鏇樺灩閻g兘顢涘☉姗嗗殼闁诲孩绋掗敋濞存粠鍨跺娲川婵犲嫮鐣垫繝娈垮灥妞存悂骞嗛弮鍫濐潊闁挎稑瀚倴濠碉紕鍋戦崐鏍礉濡ゅ懎绐楅幖娣灮椤╂彃螖閿濆懎鏆為柣鎾寸洴閺屾盯濡烽敐鍛瀴闂佹眹鍔嶉崹鍧楀蓟閿濆鍋勯柛娆忣槹閻濇棃姊虹€圭姵顥夋い锔炬暬閻涱喖螣閼测晝顦╅梺缁樏畷顒勵敆閵忊€茬箚闁绘劦浜滈埀顒佺墪鐓ゆ繝闈涙閺嬪秹鏌¢崶鈺佷憾缂傚倹宀搁悡顐﹀炊閵娧€妲堥悗鐟版啞缁诲啴濡甸崟顖氱婵°倐鍋撻柛鐕佸灦椤㈡瑩鏁撻敓锟�20濠电姴鐥夐弶搴撳亾濡や焦鍙忛柣鎴f绾惧鏌eΟ娆惧殭缂佺姴鐏氶妵鍕疀閹炬惌妫″銈庡亝濞叉ḿ鎹㈠┑瀣棃婵炴垵宕崜鎵磽娴e搫校闁搞劌娼″濠氬Χ閸℃ê寮块梺褰掑亰閸忔﹢宕戦幘婢勬棃鍩€椤掑嫬鐓濋柡鍐ㄧ墕椤懘鏌eΟ鐑橆棤闁硅櫕鎹囬妶顏呭閺夋垹顦ㄩ梺鍐叉惈閿曘儵鏁嶉崨顖滅=闁稿本鐟чˇ锔姐亜閿旇鐏︽い銏″哺椤㈡﹢濮€閻橀潧濮︽俊鐐€栧濠氬磻閹惧绡€闁逞屽墴閺屽棗顓奸崨顖ょ幢闂備胶绮濠氬储瑜斿鍛婄瑹閳ь剟寮婚弴銏犻唶婵犲灚鍔栨晥闂備胶枪妤犲摜绮旇ぐ鎺戣摕婵炴垯鍨归崡鎶芥煏婵炲灝鍔氭い顐熸櫊濮婄儤瀵煎▎鎴犳殸缂傚倸绉撮敃顏堢嵁閸愩剮鏃堝礃閳轰焦鐎梻浣告啞濞诧箓宕f惔銊ユ辈闁跨喓濮甸埛鎴︽煕濠靛棗顏い銉﹀灴閺屾稓鈧綆鍋呭畷灞炬叏婵犲啯銇濈€规洦鍋婂畷鐔煎垂椤愬诞鍥ㄢ拺闁告稑锕ラ埛鎰版煟濡ゅ啫鈻堟鐐插暣閺佹捇鎮╅搹顐g彨闂備礁鎲″ú锕傚礈濞嗘挻鍋熷ù鐓庣摠閳锋垿姊婚崼鐔恒€掔紒鐘冲哺閺屾盯骞樼€靛摜鐤勯梺璇″枓閳ь剚鏋奸弸搴ㄦ煙闁箑鏋ゆい鏃€娲樼换婵嬪閿濆棛銆愬銈嗗灥濡稓鍒掗崼銉ョ劦妞ゆ帒瀚崐鍨箾閸繄浠㈡繛鍛Ч閺岋繝鍩€椤掑嫬纭€闁绘垵妫楀▓顐︽⒑閸涘﹥澶勯柛瀣浮瀹曘儳鈧綆鍠楅悡鏇㈡煛閸ャ儱濡兼鐐瓷戞穱濠囧矗婢跺﹦浼屽┑顔硷攻濡炶棄鐣烽锕€绀嬫い鎰枎娴滄儳霉閻樺樊鍎滅紓宥嗙墪椤法鎹勯悜妯绘嫳闂佺ǹ绻戠划鎾诲蓟濞戙埄鏁冮柨婵嗘椤︺劑姊洪崫鍕闁告挾鍠栭獮鍐潨閳ь剟骞冨▎鎾搭棃婵炴垶顨呴ˉ姘辩磽閸屾瑨鍏屽┑顔炬暩閺侇噣鍨鹃幇浣圭稁婵犵數濮甸懝楣冩倷婵犲洦鐓ユ繝闈涙閸gǹ顭跨憴鍕婵﹥妞介幊锟犲Χ閸涱喚鈧儳鈹戦悙鎻掔骇闁搞劌娼¢獮濠偽旈崘鈺佺/闁荤偞绋堥崜婵嬫倶娓氣偓濮婅櫣娑甸崨顔兼锭闂傚倸瀚€氭澘鐣烽弴銏犵闁挎棁妫勯埀顒傛暬閺屻劌鈹戦崱娑扁偓妤侇殽閻愮榿缂氱紒杈ㄥ浮閹晛鐣烽崶褉鎷伴梻浣告惈婢跺洭宕滃┑鍡╁殫闁告洦鍋€濡插牊绻涢崱妤佺濞寸》鎷�
相关话题/matrix method black