# Perturbative Approach to Analog Hawking Radiation
in dielectric media in subcritical regime

S. Trevisan3,4, F. Belgiorno,1,2 S.L. Cacciatori3,4 1 Dipartimento di
Matematica, Politecnico di Milano, Piazza Leonardo 32, IT-20133 Milano, Italy
2 INdAM-GNFM  3 Department of Science and High Technology, Università
dell’Insubria, Via Valleggio 11, IT-22100 Como, Italy  4 INFN sezione di
Milano, via Celoria 16, IT-20133 Milano, Italy

###### Abstract.

We take into account the subcritical case for dielectric media by exploiting
an approximation allowing us to perform perturbative analytical calculations
and still not implying low dispersive effects. We show that in the background
of a specific soliton-like solution, pair-creation occurs and can display a
thermal behaviour governed by an effective temperature. The robustness of the
approach is also corroborated by the analysis of the $\phi\psi$-model related
to the standard Hopfield model, for which analogous results are obtained.

##  1\. Introduction

We study the problem of the so-called subcritical case for the Hawking effect
in Analogue Gravity, which is characterized by the absence of a real horizon
but still a process of pair-creation takes place. The phenomenon has been
extensively studied in the physical literature, because of the Vancouver
experiment with water [1, 2], where a thermal spectrum is detected, despite
the actual absence of a horizon, as shown in subsequent analysis [3, 4, 5]. As
to analytical calculations, the first one, again devoted to the case of
scattering in water, appeared in [6], and it was provided by means of the so-
called Bremmer approximation [7, 8], adapted to an ordinary differential
equation of the $n$-th order. Subsequent analytical calculations, involving
also the subcritical case for water and BEC, again in the limit of low
dispersive effects, are found in [9, 10], and exploit a KdV equation emerging
in the approximation of no co-propagating modes [9].
For the case of the analogue Hawking effect in dielectric media, horizonless
situations were numerically taken into account in [11, 12], but for the
subcritical case, a fully analytical calculation is still missing.
We have shown in previous papers [13, 14] that an equation of the Orr-
Sommerfeld type inherited by Nishimoto’s works (see e.g. [15, 16]) enables us
to treat in a quite unified way the case with a horizon, also called the
transcritical case, which amounts to the presence of a real turning point
(real horizon) in the limit of weak dispersive effects, where the weakness of
dispersion is indicated by a suitable small parameter $\epsilon$.
In this paper, we focus on the subcritical case for dielectric media, i.e. we
take into account configurations where a real turning point is missing. We
still obtain a fourth-order differential equation governing the phenomenon,
but we adopt a different attitude and a different expansion parameter with
respect to the transcritical case mentioned above. Indeed, we consider a
linearization of the equation around a specific soliton-like background
solution, and we exploit an expansion in terms of a parameter $\eta$ which
represents the weakness of the soliton amplitude. In the comoving frame of the
background solution, we obtain a static situation that mimics that of the
Hawking effect in the transcritical case, but with no horizon. Even if it
could seem that weak dispersion is implied by our picture, actually this is
not the case: Indeed, the dispersion parameter $\epsilon$ is a priori not
restricted to be small, and this is the main difference with respect to the
standard picture described in [13], where $\epsilon$ is small and a so-called
singular perturbation theory is to be allowed. In the framework we discuss,
one is allowed to adopt a regular perturbative expansion in the parameter
$\eta$, which represents a strong advantage.
We show that for the particular profile which is taken into account, in some
limit thermality is simulated, with an effective temperature which is one-
third of the one of the corresponding transcritical case.
Explicit calculations are firstly carried out for a modified $\phi\psi$-model
whose aim is to simplify as possible the dispersion relation associated with
the model, and to allow more straightforward analytical calculations and a
more clear exposition of the basic idea. Then we corroborate the robustness of
the approach by applying it to the $\phi\psi$ model introduced in [17] and
then discussed elsewhere and again taken into account in [13, 18]. The
aforementioned $\phi\psi$ model will be mentioned as the ‘original’ model and
represents a simplification of the Hopfield model [19] which is a standard way
to discuss the electromagnetic field in dielectric media.

##  2\. The modified $\phi\psi$-model and its solitonic solutions

The model we will consider is defined by a Lagrangian function involving two
real scalar fields $\phi$ and $\psi$, with the aim of simulating some features
of the behavior of the electromagnetic field in dielectric media. As well
known, a way to obtain this goal consists in the Hopfield model [19], where
the interaction between the electromagnetic field and the atoms/molecules of
the dielectric medium is taken into account by replacing the aforementioned
microscopic objects with a mesoscopic polarization field. One still obtains
the correct Sellmeier equation for the dispersion relation. In our case,
intending to simplify as most as possible both the analytical calculations and
the dispersion relation, we replace the polarization field with the field
$\psi$ and the electromagnetic field with $\phi$, with a set-up aimed at
reproducing the Cauchy dispersion relation most straightforwardly. We stress
that our picture below will be corroborated also in the trickier case of the
$\phi\psi$-model discussed in [17], where the Hopfield model is reduced in the
most direct way to a model reproducing exactly the Sellmeier dispersion
relation using a couple of scalar fields $\phi,\psi$.
In our present model, the Lagrangian, expressed in the lab frame, with respect
to spacetime variables $t_{l},x_{l}$, is

$\displaystyle\mathcal{L}=\frac{1}{2}(\partial_{t_{l}}\phi)^{2}+\frac{1}{2}\left((\partial_{t_{l}}\psi)^{2}+\mu^{2}\psi^{2}\right)+g\phi\partial_{x_{l}}\psi-\frac{\lambda}{4!}\psi^{4}\,.$
(2.1)

We can write it in a covariant form, which will be useful to pass to the frame
comoving with the pulse:

$\displaystyle\mathcal{L}=\frac{1}{2}(v^{\nu}\partial_{\nu}\phi)^{2}+\frac{1}{2}\left((v^{\nu}\partial_{\nu}\psi)^{2}+\mu^{2}\psi^{2}\right)+g\phi
n^{\nu}\partial_{\nu}\psi-\frac{\lambda}{4!}\psi^{4}\,,$ (2.2)

where in the lab frame we have $v^{\nu}_{lab}=(1,0)$ and
$n^{\nu}_{lab}=(0,1)$, while in a boosted frame we get
$v^{\nu}=(\gamma,-V\gamma)$ and $n^{\nu}=(-V\gamma,\gamma)$ ($c=1$
everywhere).

The equations of motion that follow from (2.2) are

$\displaystyle(v^{\nu}\partial_{\nu})^{2}\phi-gn^{\nu}\partial_{\nu}\psi=0,$
(2.3)
$\displaystyle(v^{\nu}\partial_{\nu})^{2}\psi+gn^{\nu}\partial_{\nu}\phi-\mu^{2}\phi+\frac{\lambda}{3!}\psi^{3}=0.$
(2.4)

The free-field solutions (for $\lambda=0$) are plane waves
$e^{i\omega_{lab}t_{l}-ik_{lab}x_{l}}$ which satisfy

$\displaystyle
n_{0}^{2}(\omega_{lab}):=\frac{k_{lab}^{2}}{\omega_{lab}^{2}}=\frac{\mu^{2}}{g^{2}}+\frac{\omega_{lab}^{2}}{g^{2}}=:A+B\omega_{lab}^{2}\,,$
(2.5)

where $\omega_{lab}=v^{\mu}k_{\mu}$ and $k_{lab}=-n^{\mu}k_{\mu}$. The
quantity $n_{0}$ is, by definition, the refractive index of the medium. Eq.
(2.5), in the limit as $\omega_{lab}\to 0$, gives rise to the Cauchy
dispersion relation, which is frequently used to describe dielectric media at
low frequencies. Indeed, as $\omega_{lab}\to 0$ we get

$\displaystyle
n_{0}(\omega_{lab})\simeq\sqrt{A}+\frac{B}{2\sqrt{A}}\omega_{lab}^{2}\,.$
(2.6)

To describe ordinary dielectrics, we require $A>1$, i.e. $\mu>g$.

The Lagrangian (2.2) admits a conserved current, which is related by the
Noether theorem to the invariance (in its complexified version) under phase
shifting of the fields $\phi\mapsto e^{i\alpha}\phi$, $\psi\mapsto
e^{i\alpha}\psi$:

$\displaystyle
J^{\nu}:=\frac{i}{2}\left[(v^{\nu}\phi^{\ast}(v^{\alpha}\partial_{\alpha}\phi)+v^{\nu}\psi^{\ast}(v^{\alpha}\partial_{\alpha}\psi)+gn^{\nu}\psi^{\ast}\phi-c.c.\right].$
(2.7)

The zeroth component gives rise to the conserved charge

$\displaystyle\int_{\Sigma_{t}}dxJ^{0}=:Q,$ (2.8)

where $\Sigma_{t}$ is a spacelike hypersurface. The associated (conserved)
scalar product is

$\displaystyle\left(\begin{pmatrix}\phi\\\
\psi\end{pmatrix},\begin{pmatrix}\tilde{\phi}\\\
\tilde{\psi}\end{pmatrix}\right)=\frac{i}{2}\int_{\Sigma_{t}}\left[v^{0}\phi^{\ast}v^{\alpha}\partial_{\alpha}\tilde{\phi}-v^{0}\tilde{\phi}v^{\alpha}\partial_{\alpha}\phi^{\ast}+v^{0}\psi^{\ast}v^{\alpha}\partial_{\alpha}\tilde{\psi}-v^{0}\tilde{\psi}v^{\alpha}\partial_{\alpha}\psi^{\ast}+gn^{0}(\psi^{\ast}\tilde{\phi}-\tilde{\psi}\phi^{\ast})\right].$
(2.9)

It is straightforward now to provide a canonical quantization of the free-
lagrangian, compute the propagators and perform standard perturbative QFT
computations with the vertex $i\lambda$: we will not deal with this in this
paper. The free particles of the theory are polaritons which satisfy the
Cauchy dispersion relation and the interaction term is a modelization of the
nonlinear response that can happen in dielectrics, better known as the Kerr
effect. The norm of the free particle states
$\Phi_{k}=\begin{pmatrix}\phi_{k}\\\ \psi_{k}\end{pmatrix}\sim
e^{-ik_{\mu}x^{\mu}}$ has a simple expression in the laboratory frame:

$\displaystyle||\Phi_{k}||^{2}=Q(\phi_{k},\psi_{k})\propto
2\omega_{lab}\left(1+\frac{\omega_{lab}^{2}}{\mu^{2}+\omega_{lab}^{2}}\right)\,.$
(2.10)

From Eq. (2.10) we can clearly distinguish the positive-norm modes
($\omega_{lab}>0$) from the negative-norm modes ($\omega_{lab}<0$). In
Appendix A we give a more detailed analysis of the current two-vector
$J^{\mu}$ and we prove that the norm is proportional to
$\text{sign}(\omega_{lab})$ in _any_ inertial frame.

The equations of motion (2.3) also admit a solitonic solution that propagates
rigidly at a fixed velocity $V$ with respect to the laboratory:

$\displaystyle\psi_{s}(x_{l}-Vt_{l})$
$\displaystyle=\frac{\alpha}{\cosh(\beta(x_{l}-Vt_{l}))}\,,$ (2.11)
$\displaystyle\alpha^{2}$ $\displaystyle=\frac{12V^{2}\beta^{2}}{\lambda}\,,$
(2.12) $\displaystyle\beta^{2}$
$\displaystyle=\frac{1}{V^{4}}(\mu^{2}V^{2}-g^{2})\,.$ (2.13)

##  3\. The Linearized EOM

In this section we consider the linearization of the EOMs (2.3) around a
soliton-like background:

$\displaystyle\psi_{B}(x_{l}-Vt_{l})=\sqrt{\frac{2|\eta|}{|\lambda|}}\frac{1}{\cosh(\beta(x_{l}-Vt_{l}))}\,.$
(3.1)

We will treat $\eta$ and $\beta$ as independent parameters to allow valid
results not only for the soliton (2.11) but also for other backgrounds with
the same shape. Indeed, it is known that solitons of this type can emerge also
as _approximated_ solutions of electrodynamics inside nonlinear dielectrics:
specifically, they are solutions of the so-called nonlinear Schrödinger
equation, which derives from some approximations made on the Maxwell equations
inside those media. We shall comment further on our choice in the following
section.

We want to study the scattering of the asymptotic normal modes against the
perturbation given by the soliton travelling across the medium: if the
scattering involves both positive and negative-norm modes, this can be
interpreted as a sign of instability of the system, which decays by emitting
particles, as in the pioneering computation by Hawking for black hole
evaporation.

The linearized EOMs are

$\displaystyle(v^{\nu}\partial_{\nu})^{2}\phi-gn^{\nu}\partial_{\nu}\psi=0,$
(3.2)
$\displaystyle(v^{\nu}\partial_{\nu})^{2}\psi+gn^{\nu}\partial_{\nu}\psi-\mu^{2}\phi+\frac{\lambda}{2}\psi_{B}^{2}\psi=0.$
(3.3)

Notice that, in the linearized equation, the effect of the perturbation is a
shift of the parameter $\mu^{2}$:

$\displaystyle\mu^{2}\mapsto\mu^{2}-\frac{\lambda}{2!}\psi^{2}_{B}(x_{l}-Vt_{l})\,,$
(3.4)

which induces a shift in the refractive index (see (2.5))

$\displaystyle n^{2}(\omega_{lab},x_{l},t_{l})$
$\displaystyle=n_{0}^{2}(\omega_{lab})+\delta n^{2}(x_{l}-Vt_{l})\,,$ (3.5)
$\displaystyle\delta n^{2}(x_{l}-Vt_{l})$
$\displaystyle=-\frac{\eta}{g^{2}}\frac{1}{\cosh^{2}(\beta(x_{l}-Vt_{l}))}\,.$
(3.6)

In order to include both cases $\delta n>0$ and $\delta n<0$, we will consider
$\eta$ as a real quantity, with

$\text{sign}(\eta):=\text{sign}(\lambda)\,.$

The dispersion relations (2.5) and (3.5) are represented in Figure 1, as seen
in the comoving frame with the background.

(a)

(b)

Figure 1.  (A): The dispersion relation (2.5) represented in the comoving
frame with the background, with $g=1$, $\mu=1.2$; for $0<\omega<\omega_{MAX}$
there are four real solutions. (B): The dispersion relation (3.5) represented
at the peak of the perturbation, _i.e._ $x_{l}=Vt_{l}$, for $\eta=-1$. The
modes $0<\omega<\omega_{MIN}$ do not experience an event horizon (_subcritical
regime_).

We decouple the above equations (3.2) and (3.3), by applying the operator
$(v^{\nu}\partial_{\nu})^{2}$ to the second equation; then, we eliminate
$\phi$ using $(v^{\nu}\partial_{\nu})^{2}\phi=gn^{\nu}\partial_{\nu}\psi$. In
this way, we obtain the following decoupled equation only for the field
$\psi$:

$\displaystyle(v^{\nu}\partial_{\nu})^{4}\psi+g^{2}(n^{\nu}\partial_{\nu})^{2}\psi-\mu^{2}(v^{\nu}\partial_{\nu})^{2}\phi+\frac{\lambda}{2}(v^{\nu}\partial_{\nu})^{2}(\psi_{B}^{2}\psi)=0.$
(3.7)

It is convenient to write the equations in the comoving coordinates
$t=\gamma(t_{l}-Vx_{l})$, $x=\gamma(x_{l}-Vt_{l})$. Since the potential term
is independent of the comoving time, we seek a solution in the form
$\psi=e^{-i\omega t}f(x)$, where the constant $\omega$ is a conserved quantity
(cf. e.g. [17]). In this way we end up with a fourth-order equation for $f(x)$
only:

$\displaystyle 0=$ $\displaystyle
V^{4}\gamma^{4}f^{(4)}(x)+4iV^{3}\omega\gamma^{4}f^{(3)}(x)+\gamma^{2}f^{\prime\prime}(x)\left(\eta
V^{2}\text{sech}^{2}(\beta
x/\gamma)+g^{2}-\mu^{2}V^{2}-6V^{2}\omega^{2}\gamma^{2}\right)$ (3.8)
$\displaystyle+2iV\gamma f^{\prime}(x)\left(\eta\,\text{sech}^{2}(\beta
x/\gamma)(2i\beta V\tanh(\beta
x/\gamma)+\omega\gamma)+\omega\gamma\left(g^{2}-\mu^{2}-2\omega^{2}\gamma^{2}\right)\right)$
(3.9) $\displaystyle+f(x)\left(-2\beta^{2}\eta V^{2}\text{sech}^{4}(\beta
x/\gamma)-\eta\,\text{sech}^{2}(\beta x/\gamma)(\omega\gamma+2i\beta
V\tanh(\beta
x/\gamma))^{2}+\omega^{2}\gamma^{2}\left(-g^{2}V^{2}+\mu^{2}+\omega^{2}\gamma^{2}\right)\right)\,.$
(3.10)

We can further manipulate the equation by performing the following change of
variables, usual for the Pöschl-Teller potential,

$\displaystyle z$
$\displaystyle=-e^{2\tilde{\beta}x}\,,\quad\tilde{\beta}:=\frac{\beta}{\gamma}\,,$
(3.12)

which implies

$\displaystyle\partial_{x}$
$\displaystyle=2\tilde{\beta}\theta_{z}:=2\tilde{\beta}z\frac{d}{dz}\,.$
(3.13)

By defining the rescaled parameters $G=\frac{g}{2\tilde{\beta}}$,
$\Omega=\frac{\omega}{2\tilde{\beta}}$, $M=\frac{\mu}{2\tilde{\beta}}$, we end
up with the following equation

$\displaystyle 0=$
$\displaystyle\Big{(}V^{4}\gamma^{4}\Big{)}z^{4}f^{(4)}+\Big{(}6V^{4}\gamma^{4}-4iV^{3}\Omega\gamma^{4}\Big{)}z^{3}f^{(3)}$
(3.14)
$\displaystyle+\left(G^{2}\gamma^{2}-M^{2}V^{2}\gamma^{2}+7V^{4}\gamma^{4}-12iV^{3}\Omega\gamma^{4}-6V^{2}\Omega^{2}\gamma^{4}-\frac{\eta}{\tilde{\beta}^{2}}\frac{V^{2}\gamma^{2}z}{(1-z)^{2}}\right)z^{2}f^{(2)}$
(3.15)
$\displaystyle+\Bigg{(}G^{2}\gamma^{2}-M^{2}V^{2}\gamma^{2}-2iG^{2}V\Omega\gamma^{2}+2iM^{2}V\Omega\gamma^{2}+V^{4}\gamma^{4}-4iV^{3}\Omega\gamma^{4}-6V^{2}\Omega^{2}\gamma^{4}+4iV\Omega^{3}\gamma^{4}$
(3.16)
$\displaystyle+\frac{\eta}{\tilde{\beta}^{2}}\left(\frac{2iV\Omega\gamma^{2}z}{(1-z)^{2}}-\frac{V^{2}\gamma^{2}z(3+z)}{(1-z)^{3}}\right)\Bigg{)}\,z\,f^{(1)}$
(3.17)
$\displaystyle+\left(M^{2}\Omega^{2}\gamma^{2}-G^{2}V^{2}\Omega^{2}\gamma^{2}+\Omega^{4}\gamma^{4}+\frac{\eta}{\tilde{\beta}^{2}}\left(\frac{\Omega^{2}\gamma^{2}z}{(1-z)^{2}}+\frac{2iV\Omega\gamma^{2}z}{(1-z)^{3}}-\frac{V^{2}\gamma^{2}z(1+4z+z^{2})}{(1-z)^{4}}\right)\right)f\,.$
(3.18)

Equation (3.14) is of Fuchsian type. The complete characterization of the
equation is not in the interest of this work and will be the subject of some
future publication. In the present work, we propose a perturbative method for
the solution of the equation, which will be exposed in the next Section.

##  4\. Perturbative Method

The solitonic solution (2.11) has no free parameter, except for the velocity
$V$: its amplitude, in particular, is fixed by the constants of the
Lagrangian. Nevertheless, as in the previous section, we linearize the field
equations for our nonlinear theory around weak soliton-like backgrounds
(whichever produced experimentally), with the aim to consider the perturbation
induced by the soliton-like background as “small”. We stress that we have the
freedom to choose the spatial dependence of the perturbation in the refractive
index in such a way that it can be reproduced using linearization around a
given background solution. This attitude is corroborated by the experimental
fact that strong laser pulses may induce in dielectrics a non-linear
perturbation $\delta n(x)$ (in the comoving frame) of the refractive index
which is typically orders of magnitude lower than the leading term (e.g. order
of $10^{-3}$ compared with one of the leading term).
There is an immediate consequence of this approach: as all the non-homogeneity
is associated with the refractive index $n(x)$ is due to the correction
$\delta n(x)$, and the latter one is treated perturbatively in $\eta$ which is
associated directly with the amplitude of the soliton-like background, we are
considering a situation where at the lowest order homogeneity occurs, and then
no horizon can appear, i.e. we are automatically in a subcritical regime. This
is the main difference with respect to the perturbative approach where the
perturbation parameter $\epsilon$ is associated with a weak dispersion, and
where instead a horizon may appear in the leading order equation (in form of a
real turning point).

The perturbative method for Fuchsian equations that we apply was proposed in
[20]. We start by eliminating the third order derivative by the change of
variable

$f(z)=z^{-\frac{3}{2}+i\frac{\Omega}{V}}u(z)\,,$

that gives

$\displaystyle u^{(4)}+(v_{1}(z)+\eta w_{1}(z))u^{(2)}+(v_{2}(z)+\eta
w_{2}(z))u^{(1)}+(v_{3}(z)+\eta w_{3}(z))u=0\,,$ (4.1)

where

$\displaystyle v_{1}(z)$
$\displaystyle=\frac{2G^{2}-2M^{2}V^{2}+5V^{4}\gamma^{2}}{2V^{4}\gamma^{2}z^{2}}\,,$
$\displaystyle v_{2}(z)$
$\displaystyle=\frac{-2G^{2}V+2M^{2}V^{3}+2iG^{2}\Omega-2iG^{2}V^{2}\Omega-5V^{5}\gamma^{2}}{V^{5}\gamma^{2}z^{3}}\,,$
$\displaystyle v_{3}(z)$
$\displaystyle=\frac{36G^{2}V^{2}-36M^{2}V^{4}-48iG^{2}V\Omega+48iG^{2}V^{3}\Omega-16G^{2}\Omega^{2}+32G^{2}V^{2}\Omega^{2}-16G^{2}V^{4}\Omega^{2}+81V^{6}\gamma^{2}}{16V^{6}\gamma^{2}z^{4}}\,,$
$\displaystyle w_{1}(z)$
$\displaystyle=-\frac{1}{\tilde{\beta}^{2}V^{2}\gamma^{2}(z-1)^{2}z}\,,$
$\displaystyle w_{2}(z)$
$\displaystyle=\frac{4}{\tilde{\beta}^{2}V^{2}\gamma^{2}(z-1)^{3}z}\,,$
$\displaystyle w_{3}(z)$
$\displaystyle=\frac{-1+2z-25z^{2}}{4\tilde{\beta}^{2}V^{2}\gamma^{2}(z-1)^{4}z^{3}}\,.$

Now we consider $\eta$ as a small parameter, and formally expand the solution
as

$\displaystyle u(z)=u_{0}(z)+\eta u_{1}(z)+\eta^{2}u_{2}(z)+...$ (4.2)

which allows us to obtain a regular perturbative expansion (to be compared
with the singular perturbation expansion one obtains by expanding with respect
to a low dispersion parameter $\epsilon$, cf. section 7). Herein, we compute
the first order solution $u_{1}$ and we will discuss the possibility of taking
$\eta\rightarrow 1$.

Solving the unperturbed equation ($\eta=0$) is very easy and it gives

$\displaystyle u_{0}(z)=z^{i\alpha}\,,$ (4.3)

where $\alpha$ satisfies a fourth-degree algebraic equation. By putting
$i\alpha=i\frac{k}{2\beta}+\frac{3}{2}+i\frac{\Omega}{V}$, we find that $k$ is
one of the four solutions of the dispersion relation (2.5) as written in the
comoving frame:

$\displaystyle\gamma\left(-g^{2}(k+V\omega)^{2}+(kV+\omega)^{2}(\mu^{2}+(kV+\omega)^{2}\gamma^{2})\right)=0\,.$
(4.4)

By substituting (4.2) into (4.1) we find the set of equations

$\displaystyle
u_{1}^{(4)}+v_{1}(z)u_{1}^{(2)}+v_{2}(z)u_{1}^{(1)}+v_{3}(z)u_{1}$
$\displaystyle=w_{1}(z)u_{0}^{(2)}+w_{2}(z)u_{0}^{(1)}+w_{3}(z)u_{0}\,,$ (4.5)
$\displaystyle...\,,$ (4.6) $\displaystyle
u_{n}^{(4)}+v_{1}(z)u_{n}^{(2)}+v_{2}(z)u_{n}^{(1)}+v_{3}(z)u_{n}$
$\displaystyle=w_{1}(z)u_{(n-1)}^{(2)}+w_{2}(z)u_{(n-1)}^{(1)}+w_{3}(z)u_{(n-1)}$
(4.7)

Thus every $u_{n}$ satisfies a linear differential equation with a source term
that depends on $u_{(n-1)}$; the associated homogeneous equation is the
unperturbed equation satisfied by $u_{0}$.

Let us consider the equation for $u_{1}$. We can explicitly solve it by
applying the method of variation of constants. By defining

$\displaystyle
r_{1}(z):=w_{1}(z)u_{0}^{(2)}+w_{2}(z)u_{0}^{(1)}+w_{3}(z)u_{0}\,,$ (4.8)

we compute the quantities

$\displaystyle W(z)$
$\displaystyle=\det\begin{pmatrix}z^{i\alpha_{1}}&z^{i\alpha_{2}}&z^{i\alpha_{3}}&z^{i\alpha_{4}}\\\
(z^{i\alpha_{1}})^{\prime}&(z^{i\alpha_{2}})^{\prime}&(z^{i\alpha_{3}})^{\prime}&(z^{i\alpha_{4}})^{\prime}\\\
(z^{i\alpha_{1}})^{\prime\prime}&(z^{i\alpha_{2}})^{\prime\prime}&(z^{i\alpha_{3}})^{\prime\prime}&(z^{i\alpha_{4}})^{\prime\prime}\\\
(z^{i\alpha_{1}})^{\prime\prime\prime}&(z^{i\alpha_{2}})^{\prime\prime\prime}&(z^{i\alpha_{3}})^{\prime\prime\prime}&(z^{i\alpha_{4}})^{\prime\prime\prime}\\\
\end{pmatrix}\,,$ $\displaystyle W_{1}(z)$
$\displaystyle=r_{1}(z)\det\begin{pmatrix}z^{i\alpha_{2}}&z^{i\alpha_{3}}&z^{i\alpha_{4}}\\\
(z^{i\alpha_{2}})^{\prime}&(z^{i\alpha_{3}})^{\prime}&(z^{i\alpha_{4}})^{\prime}\\\
(z^{i\alpha_{2}})^{\prime\prime}&(z^{i\alpha_{3}})^{\prime\prime}&(z^{i\alpha_{4}})^{\prime\prime}\\\
\end{pmatrix}\,,$ $\displaystyle W_{2}(z)$
$\displaystyle=-r_{1}(z)\det\begin{pmatrix}z^{i\alpha_{1}}&z^{i\alpha_{3}}&z^{i\alpha_{4}}\\\
(z^{i\alpha_{1}})^{\prime}&(z^{i\alpha_{3}})^{\prime}&(z^{i\alpha_{4}})^{\prime}\\\
(z^{i\alpha_{1}})^{\prime\prime}&(z^{i\alpha_{3}})^{\prime\prime}&(z^{i\alpha_{4}})^{\prime\prime}\\\
\end{pmatrix}\,,$ $\displaystyle W_{3}(z)$
$\displaystyle=r_{1}(z)\det\begin{pmatrix}z^{i\alpha_{1}}&z^{i\alpha_{2}}&z^{i\alpha_{4}}\\\
(z^{i\alpha_{1}})^{\prime}&(z^{i\alpha_{2}})^{\prime}&(z^{i\alpha_{4}})^{\prime}\\\
(z^{i\alpha_{1}})^{\prime\prime}&(z^{i\alpha_{2}})^{\prime\prime}&(z^{i\alpha_{4}})^{\prime\prime}\\\
\end{pmatrix}\,,$ $\displaystyle W_{4}(z)$
$\displaystyle=-r_{1}(z)\det\begin{pmatrix}z^{i\alpha_{1}}&z^{i\alpha_{2}}&z^{i\alpha_{3}}\\\
(z^{i\alpha_{1}})^{\prime}&(z^{i\alpha_{2}})^{\prime}&(z^{i\alpha_{3}})^{\prime}\\\
(z^{i\alpha_{1}})^{\prime\prime}&(z^{i\alpha_{2}})^{\prime\prime}&(z^{i\alpha_{3}})^{\prime\prime}\\\
\end{pmatrix}\,.$

A particular solution to the first equation of (4.5) is given by

$\displaystyle u_{1}(z)=z^{i\alpha_{1}}\int
dz^{\prime}\frac{W_{1}(z^{\prime})}{W(z^{\prime})}+z^{i\alpha_{2}}\int
dz^{\prime}\frac{W_{2}(z^{\prime})}{W(z^{\prime})}+z^{i\alpha_{3}}\int
dz^{\prime}\frac{W_{3}(z^{\prime})}{W(z^{\prime})}+z^{i\alpha_{4}}\int
dz^{\prime}\frac{W_{4}(z^{\prime})}{W(z^{\prime})}\,.$ (4.9)

In principle, we can recursively iterate the procedure to obtain particular
solutions to $u_{n}$. Notice that the general solution for each order $n$ is
obtained by adding a solution to the homogeneous equation, which is the same
for every $n$, so we can say that the general solution $u(z)$ is obtained by
adding to the iterative solution a combination of the solutions to the
unperturbed equation: the coefficients of the combination will be set by the
boundary condition of the scattering.

##  5\. Boundary Conditions for the Scattering

For the subcritical case, it is standard to consider a white hole
configuration, in which an initial state, representing a Hawking mode,
approaches the white hole-like perturbation at early times, and four modes
emerge at late times: three backward modes $P,N,B$, where $N$ is the only
negative norm mode appearing in the scattering, and a transmitted-mode $T$
representing the fraction of the Hawking mode $H$ which is transmitted beyond
the perturbation This is nothing but what happens in presence of a real white
hole horizon in the transcritical case, apart for the transmitted mode. For
the black hole case, an analytical study is trickier, as the black hole-like
configuration, in principle, is not related to the white hole-like one by time
reversal: in presence of a black hole horizon, one has three entering initial
modes $P,N,B$ which are converted in a scattered emerging mode $H$. When the
scattering is subcritical, in principle, given three initial modes $P,N,B$,
one should consider the possibility to get three transmitted particles, and
then it is evident that one does not obtain the time reversal configuration of
a white hole subcritical scattering.

We consider a white hole-like configuration, as usual. Thus, in the initial
state, we will have only a right-moving $H$-mode, which is scattered and mode
converted in a way such that we obtain the four modes in the final state. The
asymptotic form of the solution we seek is thus

$\displaystyle
f(x)\sim\begin{cases}e^{ik_{H}x}+C_{2}e^{ik_{P}x}+C_{3}e^{ik_{N}x}+C_{4}e^{ik_{B}x}\,,\quad&x\rightarrow-\infty\\\
C_{1}e^{ik_{H}x}\,,\quad&x\rightarrow+\infty\end{cases}\,,$ (5.1)

when we recall that the $H$-mode has a positive group velocity (right-moving),
while all other modes have negative group velocity (left-moving). Since the
unperturbed solution would be just $f(x)=e^{ik_{H}x}$ defined everywhere, we
expect the coefficient $C_{1}$ of the transmitted part of the Hawking mode to
be $O(1)$, while all other coefficients must be $O(\varepsilon)$. We can now
proceed to the solution of the equation (4.5) for $u_{1}(z)$ considering for
the source term the unperturbed solution $u_{0}(z)=z^{i\alpha_{H}}$, where
$i\alpha_{H}=i\frac{k_{H}}{2\beta}+\frac{3}{2}+i\frac{\Omega}{V}$. The
detailed computation of the solution is shown in Appendix C: in what follows
we will just expose the final expressions.

The asymptotic coefficients of the solution (see Eq. (5.1)) at first
perturbative order are

$\displaystyle C_{1}$ $\displaystyle=\frac{1}{1+i\eta\tau}\,,$ (5.2)
$\displaystyle C_{2}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\,\frac{\pi(k_{P}V+\omega)^{2}(-1)^{i(k_{H}-k_{P})}}{4\beta^{2}V^{4}(k_{P}-k_{N})(k_{P}-k_{B})}\mathrm{csch}\left((k_{H}-k_{P})\frac{\pi\gamma}{2\beta}\right)\,,$
(5.3) $\displaystyle C_{3}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\,\frac{\pi(k_{N}V+\omega)^{2}(-1)^{i(k_{H}-k_{N})}}{4\beta^{2}V^{4}(k_{N}-k_{P})(k_{N}-k_{B})}\mathrm{csch}\left((k_{H}-k_{N})\frac{\pi\gamma}{2\beta}\right)\,,$
(5.4) $\displaystyle C_{4}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\,\frac{\pi(k_{B}V+\omega)^{2}(-1)^{i(k_{H}-k_{B})}}{4\beta^{2}V^{4}(k_{B}-k_{P})(k_{B}-k_{N})}\mathrm{csch}\left((k_{H}-k_{B})\frac{\pi\gamma}{2\beta}\right)\,,$
(5.5) $\displaystyle\tau$ $\displaystyle=\frac{(k_{H}V+\omega)^{2}}{2\beta
V^{4}\gamma(k_{H}-k_{P})(k_{H}-k_{N})(k_{H}-k_{B})}\,.$ (5.6)

##  6\. Results

Current conservation implies that the following equation holds true:

$|J_{H}|=|J_{T}|+|J_{P}|-|J_{N}|+|J_{B}|,$ (6.1)

so that, by defining

$\displaystyle|T|$ $\displaystyle=\frac{|J_{T}|}{|J_{H}|},$ (6.2)
$\displaystyle|N|$ $\displaystyle=\frac{|J_{N}|}{|J_{H}|},$ (6.3)
$\displaystyle|B|$ $\displaystyle=\frac{|J_{B}|}{|J_{H}|},$ (6.4)

we get

$\displaystyle 1=|T|+|P|-|N|+|B|\,,$ (6.5)

where for the model labeled by $K$ we have
$|K|:=|C_{K}|^{2}\left|\frac{v_{K}\partial_{\omega}DR|_{k_{K}}}{v_{H}\partial_{\omega}DR|_{k_{H}}}\right|$,
$v_{K}(\omega)$ is the group velocity of the $K$-mode and
$\text{DR}(\omega,k)$ is the dispersion relation function defined in Eq.
(A.3): we provide a derivation of (6.5) in Appendix B. The ratio $|N|$
corresponds, in the black hole case, to the rate of spontaneous emission of
$H$-mode waves, as it can be argued by using standard Bogoliubov coefficients.
In the present case, it represents the rate of pair-production in the
subcritical process at hand.
From the scattering coefficients (5.2) and the expression (A.10) for the flux
factors, we can write down explicitly the expressions

$\displaystyle|N|$
$\displaystyle=\frac{\eta^{2}}{1+\eta^{2}\tau^{2}}\frac{\pi^{2}(k_{N}V+\omega)^{2}(k_{H}V+\omega)^{2}\mathrm{csch}[((k_{H}-k_{N})\gamma\pi)/(2\beta)]^{2}}{16\beta^{4}V^{8}(k_{N}-k_{P})(k_{N}-k_{B})(k_{H}-k_{P})(k_{H}-k_{B})}$
(6.6) $\displaystyle\frac{|P|}{|N|}$
$\displaystyle=\frac{(k_{N}-k_{B})(k_{P}-k_{H})(k_{P}V+\omega)^{2}(\mathrm{csch}[\pi\gamma(k_{H}-k_{P})/(2\beta)])^{2}}{(k_{P}-k_{B})(k_{N}-k_{H})(k_{N}V+\omega)^{2}(\mathrm{csch}[\pi\gamma(k_{H}-k_{N})/(2\beta)])^{2}}\,,$
(6.7) $\displaystyle|T|$ $\displaystyle=\frac{1}{1+\eta^{2}\tau^{2}}\,,$ (6.8)
$\displaystyle|B|$
$\displaystyle=\frac{\eta^{2}}{1+\eta^{2}\tau^{2}}\frac{\pi^{2}(k_{B}V+\omega)^{2}(k_{H}V+\omega)^{2}(\mathrm{csch}[\pi\gamma(k_{H}-k_{B})/(2\beta)])^{2}}{16\beta^{4}V^{8}(k_{B}-k_{N})(k_{B}-k_{P})(k_{H}-k_{N})(k_{H}-k_{P})}\,.$
(6.9)

As a consequence, we can provide

$\displaystyle|N|=\frac{1-|T|-|B|}{\frac{|P|}{|N|}-1}\,.$ (6.10)

We stress that these expressions are exact except for the perturbative
approximation in $\eta$: no other approximation has been made throughout the
computation, so they are valid for all frequencies, as long as we know
$k_{j}(\omega)$.

(a)

(b)

Figure 2. The coefficients $T$ and $B$ from (6.8) and (6.9), for different
values of the amplitude of the perturbation $\eta$.

(a)

(b)

(c)

(d)

Figure 3. The coefficients $P$ and $N$ from (6.6) and (6.7), for different
values of the amplitude of the perturbation $\eta$. In Figures (A)-(C) the
behaviour at small frequencies is $P,N\sim\omega$, which agrees with the
findings of [6]; in Figure (C) we show the onset of critical behaviour, as
$P,N\sim\omega^{-1}$.

(a)

(b)

Figure 4. The temperature $T_{\omega}$ as defined in (6.16). (A): Plot at
fixed $\beta=1$, for different $\eta$; we see the formation of a plateau at
$T_{\omega}\approx T_{pert}$ for high $\eta$. (B): Plot at fixed $\eta=1$, for
different $\beta$; the estimation of the value of the plateau $T_{pert}$ is
less accurate for increasing $\beta$.  Figure 5. The ratio
$\omega/{\log(\frac{|P|}{|N|})}$ normalized by the estimated temperature
$T_{pert}$ (see Eq. (6.15)), as a function of $\beta$. The ratio is almost
constant, which means that $\frac{|P|}{|N|}$ is indeed exponential with the
frequency; our estimation of the ”temperature” $T_{pert}$ proves to be good at
low $\beta$, while less accurate for increasing $\beta$.

For an explicit evaluation we have expanded the modes $k_{j}(\omega)$ from Eq.
(2.5) near $\omega=0$, obtaining

$\displaystyle k_{H}$ $\displaystyle=\frac{\mu-gV}{g-\mu
V}\,\omega+O(\omega^{3}),$ (6.11) $\displaystyle k_{B}$
$\displaystyle=-\frac{\mu+gV}{g+\mu V}\,\omega+O(\omega^{3}),$ (6.12)
$\displaystyle k_{P}$ $\displaystyle=\frac{\sqrt{g^{2}-\mu^{2}V^{2}}}{\gamma
V^{2}}-\left(\frac{1}{V}+\frac{g^{2}}{\gamma^{2}V(g^{2}-\mu^{2}V^{2})}\right)\,\omega-\frac{g^{2}(2g^{2}+\mu^{2}V^{2})}{2\gamma(g^{2}-\mu^{2}V^{2})^{5/2}}\,\omega^{2}+O(\omega^{3}),$
(6.13) $\displaystyle k_{N}$
$\displaystyle=-\frac{\sqrt{g^{2}-\mu^{2}V^{2}}}{\gamma
V^{2}}-\left(\frac{1}{V}+\frac{g^{2}}{\gamma^{2}V(g^{2}-\mu^{2}V^{2})}\right)\,\omega+\frac{g^{2}(2g^{2}+\mu^{2}V^{2})}{2\gamma(g^{2}-\mu^{2}V^{2})^{5/2}}\,\omega^{2}+O(\omega^{3})\,,$
(6.14)

and we put $g=1$, $\mu=1.2$, $V=0.5$. As to $k_{T}$, due to the fact that the
perturbation vanishes very rapidly in the comoving frame both for $x\to\infty$
and for $x\to-\infty$, one easily find $k_{T}=k_{H}$. These values are taken
just to obtain qualitatively the same dispersion relation as in an
experimental situation of laser pulses in silica, but they are not meant to be
quantitatively accurate. In the following analysis of the results we will
treat both $\eta$ and $\beta$ as independent parameters: although they are
uniquely determined by the soliton (2.11), we still think that it is
interesting to study the dependence of the radiation on these parameters. We
might equivalently say that we are studying a generic perturbation of the form
(3.1) with generic parameters, even though this perturbation might not be a
solution of the nonlinear equations of motion (2.3). The scattering
coefficients $|T|$, $|B|$, $|P|$ and $|N|$ are plotted in Figure 2 and Figure
3, for different values of $\eta$. Notice that $|B|$ results to be smaller
than the others, as was expected. The transmission coefficient $|T|$ is $\sim
1$ for low frequencies but it decreases more rapidly as $\eta$ increases: this
is a sign that for high perturbation we approximate the formation of an event
horizon: this interpretation is supported by the fact that the frequencies
near $\omega_{MAX}$ are less transmitted.

The coefficients $|P|$ and $|N|$ are nearly equal at low frequencies and are
$\sim\omega$. This behaviour is in agreement with what was found by Coutant
and Weinfurtner [6] for a subcritical flow in shallow water.

The ratio $\frac{|P|}{|N|}$ is very close to an exponential function,
$\frac{|P|}{|N|}\sim e^{\frac{\omega}{T}}$. This in shown if Figure 5, where
we see that $\omega/\log(\frac{|P|}{|N|})$ is almost constant. The
“temperature” $T_{pert}$ (the subscript is for “perturbative”) can be
estimated analytically from the arguments of the $\mathrm{csch}$ in (6.7): we
find

$\displaystyle\beta_{pert}=\frac{1}{T_{pert}}=\frac{2\pi\gamma}{2\beta}\lim\limits_{\omega\rightarrow
0}\frac{2k_{H}-k_{P}-k_{N}}{\omega}=\frac{2\pi g(2g+\mu V)}{\beta\gamma
V(g^{2}-\mu^{2}V^{2})}\,.$ (6.15)

We note that $T_{pert}\propto\beta$, which is proportional to the derivative
of the background function (2.11): this is what was expected based on previous
literature about analogue Hawking radiation in the critical case. As it is
clear from Figure 5, while the exponential approximation is still valid, our
estimation of the temperature is less good as $\beta$ increases.

Finally, we define a temperature function $T_{\omega}$ by the relation

$\displaystyle|N|=:\frac{1}{e^{\frac{\omega}{T_{\omega}}}-1}\,.$ (6.16)

The function $T_{\omega}$ is plotted in Figure 4 for different values of
$\eta$ and $\beta$. We note $T_{\omega}$ is not constant, showing a lack of
thermality in the emission spectrum. However, for higher values of $\eta$, a
plateau is created for frequencies close to $\omega_{MAX}$, with
$T_{\omega}\approx T_{H}$. This fact is again interpreted as a sign of the
presence of a group horizon for higher frequencies: indeed, a similar
behaviour was found by a numerical study of the transcritical regime in
shallow water [3].

The near criticality of the system for high $\eta$ is underlined also by the
behaviour of $|P|$ and $|N|$, which start growing as $\sim\frac{1}{\omega}$
near $\omega=0$ (see Figure 3(d)). This is precisely the behaviour that one
expects in the critical case. We specify that the results for high values of
$\eta$ should be taken carefully since $\eta$ is precisely the expansion
parameter of our perturbative solution. However we find it interesting that
already at the leading order, the solution shows a near critical behaviour for
high $\eta$; we expect that the same qualitative behaviour is present also in
the higher order solutions, possibly showing up already at lower values of
$\eta$.

##  7\. The original $\phi\psi$-model

In the original $\phi\psi$-model, including a $\psi^{4}$ term, one has

$\displaystyle{\mathcal{L}}_{\varphi\psi}=\frac{1}{2}(\partial_{\mu}\phi)(\partial^{\mu}\phi)+\frac{1}{2\chi\omega_{0}^{2}}\left[(v^{\alpha}\partial_{\alpha}\psi)^{2}-\omega_{0}^{2}\psi^{2}\right]+\frac{1}{c}(v^{\alpha}\partial_{\alpha}\psi)\phi-\frac{\lambda}{4!}\psi^{4},$
(7.1)

where $\phi,\psi$ play the role of the electromagnetic field and polarization
field respectively, $\chi$ plays the role of the dielectric susceptibility,
$v^{\mu}$ is the usual four-velocity vector of the dielectric, $\omega_{0}$ is
the proper frequency of the medium [18]. We get the system

$\displaystyle\Box\phi-\frac{1}{c}(v^{\mu}\partial_{\mu}\psi)$
$\displaystyle=0,$ (7.2)
$\displaystyle\left(\frac{1}{\chi\omega_{0}^{2}}(v^{\mu}\partial_{\mu})^{2}+\frac{1}{\chi}\right)\psi+\frac{1}{c}(v^{\mu}\partial_{\mu}\phi)+\frac{\lambda}{3!}\psi^{3}$
$\displaystyle=0.$ (7.3)

As in [21], we allow the spatial dependence to appear in $\chi$ and in
$\omega_{0}$ in such a way that $\chi\omega_{0}^{2}=$ const. Indeed, by
linearizing the model around soliton-like solutions, directly written in the
comoving frame, one obtains

$\displaystyle\frac{1}{\chi}$ $\displaystyle\mapsto$
$\displaystyle\frac{1}{\chi}+\frac{\lambda}{2}\psi_{0}^{2}(x),$ (7.4)
$\displaystyle\omega_{0}^{2}$ $\displaystyle\mapsto$
$\displaystyle\omega_{0}^{2}(1+\chi\frac{\lambda}{2}\psi_{0}^{2}(x)),$ (7.5)

in such a way that $\chi\omega_{0}^{2}$ remains invariant:
$\chi\omega_{0}^{2}\mapsto\chi\omega_{0}^{2}$.

In this case, we identify the parameter associated with dispersion as follows:

$\epsilon^{2}\coloneqq\frac{1}{\chi\omega_{0}^{2}}.$

This parameter has been considered as the small parameter associated with the
model in the limit of low dispersion.
In the present case, we adopt a different view where dispersion can also be
strong, and the expansion parameter is instead associated with the amplitude
of the background solution around which the EOMs are linearized. We limit
ourselves to notice that a linearization of the EOMs around a background
soliton-like solution $\psi_{0}(x)$ amounts simply in replacing
$\frac{\lambda}{3!}\psi^{3}\mapsto\frac{\lambda}{2}\psi_{0}^{2}\psi$ in (7.3).

###  7.1. A separated equation for $\psi$

By applying the operator $\Box$ on the left of equation (7.3), as shown in
[13], we obtain the following fourth-order ordinary differential equation:

$\begin{split}&-\epsilon^{2}\partial_{x}^{4}f-2i\epsilon^{2}\frac{\omega}{v}\partial_{x}^{3}f+\frac{1}{\chi\gamma^{2}v^{2}}\left(-\biggl{(}1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}\biggr{)}+\epsilon^{2}\chi\omega^{2}\right)\partial_{x}^{2}f+2\biggl{(}i\frac{\omega}{v}\frac{1}{c^{2}}(1-\epsilon^{2}\omega^{2})-\frac{1}{\gamma^{2}v^{2}}\biggl{(}\partial_{x}\frac{1}{\chi}\biggr{)}\biggr{)}\partial_{x}f\\\
&+\left(\epsilon^{2}\frac{\omega^{4}}{v^{2}c^{2}}-\frac{1}{\gamma^{2}v^{2}}\biggl{(}\partial_{x}^{2}\frac{1}{\chi}\biggr{)}-\frac{\omega^{2}}{\chi\gamma^{2}v^{2}c^{2}}-\frac{\omega^{2}}{c^{2}v^{2}}\right)f=0.\end{split}$
(7.6)

We also define $f(x)=h(x)\zeta(x)$, with

$h(x)=A\exp(-i\frac{\omega}{2v}x),$ (7.7)

where $A$ is a constant. $h(x)$ is chosen such that the third-order term
vanishes, and the procedure is analogous to the Liouville transformation which
eliminates the first-order term in a second-order linear ordinary differential
equation. This leads to the following quartic equation, which is just of the
type ‘Orr–Sommerfeld’

$\begin{split}&-\epsilon^{2}\partial_{x}^{4}\zeta+\biggl{[}-\frac{1}{\chi\gamma^{2}v^{2}}\biggl{(}1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}\biggr{)}+\epsilon^{2}\frac{1}{\gamma^{2}v^{2}}\biggl{(}1-\frac{3}{2}\gamma^{2}\biggr{)}\omega^{2}\biggr{]}\partial_{x}^{2}\zeta\\\
&+\biggl{(}i\frac{\omega}{v}\frac{1}{\chi\gamma^{2}v^{2}}\biggl{(}1+\chi\gamma^{2}\frac{v^{2}}{c^{2}}\biggr{)}-2\frac{1}{\gamma^{2}v^{2}}\biggl{(}\partial_{x}\frac{1}{\chi}\biggr{)}-i\epsilon^{2}\frac{\omega^{3}}{vc^{2}}\biggr{)}\partial_{x}\zeta\\\
&+\biggl{[}\frac{1}{\gamma^{2}v^{2}}\biggl{(}i\frac{\omega}{v}\biggl{(}\partial_{x}\frac{1}{\chi}\biggr{)}-\biggr{(}\partial_{x}^{2}\frac{1}{\chi}\biggr{)}\biggr{)}+\frac{1}{\gamma^{2}v^{2}}\biggl{(}\frac{1}{4\chi}\frac{\omega^{2}}{v^{2}}\biggl{(}1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}\biggr{)}-\frac{\omega^{2}}{\chi
c^{2}}\biggl{)}\\\
&+\epsilon^{2}\biggl{(}\frac{\omega^{4}}{v^{4}}\biggl{(}-\frac{1}{16}+\frac{1v^{2}}{4c^{2}}\biggr{)}\biggr{)}\biggr{]}\zeta=0.\end{split}$
(7.8)

The effect of the linearization around a background soliton-like solution
$\psi_{0}(x)$ consists simply in the replacement

$\displaystyle\frac{1}{\chi}\mapsto\frac{1}{\chi_{0}}+\frac{\lambda}{2}\psi_{0}^{2}(x),$
(7.9)

where $\chi_{0}$ is constant. We choose

$\displaystyle\psi_{0}(x):=2\sqrt{\frac{|\eta|}{|\lambda|}}\frac{1}{\cosh(\tilde{\beta}x)}.$
(7.10)

As in the model discussed in the previous sections, we turn to the change of
variable

$\displaystyle z:=-\exp(2\tilde{\beta}x).$ (7.11)

As a consequence, with some abuse of language, we have $\zeta=f(z)$. Then we
obtain

$\displaystyle\frac{1}{\chi}\mapsto\frac{1}{\chi_{0}}-2\eta\frac{z}{(1-z)^{2}}.$
(7.12)

Furthermore, to cancel a new third-order term that appears again after the
above-mentioned independent variable change, we set $\zeta(z)=z^{-3/2}h(z)$.
Then we obtain the following equation:

$\displaystyle h^{(4)}+(u_{2}+\eta w_{2})h^{(2)}+(u_{1}+\eta
w_{1})h^{(1)}+(u_{0}+\eta w_{0})h=0,$ (7.13)

where

$\displaystyle u_{2}$
$\displaystyle=\frac{\frac{1}{\gamma^{2}v^{2}\chi_{0}}(1-\chi_{0}\gamma^{2}\frac{v^{2}}{c^{2}})+\epsilon^{2}(10\tilde{\beta}^{2}-(1-\frac{3}{2}\gamma^{2})\omega^{2})}{4\tilde{\beta}^{2}\epsilon^{2}z^{2}}\,,$
$\displaystyle u_{1}$
$\displaystyle=\frac{-4\tilde{\beta}vc^{2}(1-\chi_{0}\gamma^{2}\frac{v^{2}}{c^{2}})-ic^{2}\omega(1+\chi_{0}\gamma^{2}\frac{v^{2}}{c^{2}})+\epsilon^{2}\chi_{0}v(-40\tilde{\beta}^{3}v^{2}c^{2}\gamma^{2}+2\tilde{\beta}c^{2}(2-3\gamma^{2})\omega^{2}+i\gamma^{2}v\omega^{3})}{8\tilde{\beta}^{3}c^{2}\chi_{0}\epsilon^{2}\gamma^{2}v^{3}z^{3}}\,,$
$\displaystyle u_{0}$
$\displaystyle=\frac{144\tilde{\beta}^{2}v^{2}c^{2}(1-\chi_{0}\gamma^{2}\frac{v^{2}}{c^{2}})+48i\tilde{\beta}vc^{2}(1+\chi_{0}\gamma^{2}\frac{v^{2}}{c^{2}})+4\omega^{2}((4+\chi_{0}\gamma^{2})v^{2}-c^{2})+\epsilon^{2}\delta}{256\tilde{\beta}^{4}c^{2}\chi_{0}\epsilon^{2}\gamma^{2}v^{4}z^{4}}\,,$
$\displaystyle w_{2}$
$\displaystyle=-\frac{1}{2\tilde{\beta}^{2}\epsilon^{2}\gamma^{2}v^{2}z(1-z)^{2}}\,,$
$\displaystyle w_{1}$
$\displaystyle=\frac{i\omega(z-1)+8\tilde{\beta}vz}{4\tilde{\beta}^{3}\epsilon^{2}\gamma^{2}v^{3}z^{2}(-1+z)^{3}}\,,$
$\displaystyle w_{0}$
$\displaystyle=\frac{(c^{2}-4v^{2})\omega^{2}(z-1)^{2}+4i\tilde{\beta}c^{2}v\omega(1-z)(5z-1)-4\tilde{\beta}^{2}c^{2}v^{2}(1-2z+25z^{2})}{32\tilde{\beta}^{4}c^{2}\epsilon^{2}\gamma^{2}v^{4}(z-1)^{4}z^{3}}\,,$

with

$\delta:=\chi_{0}(1296\tilde{\beta}^{4}v^{4}c^{2}\gamma^{2}+72\tilde{\beta}^{2}c^{2}v^{2}\omega^{2}(-2+3\gamma^{2})-48i\tilde{\beta}\gamma^{2}v^{3}\omega^{3}+\gamma^{2}(c^{2}-4v^{2})\omega^{4}).$

###  7.2. Dispersion relation and its roots

As in the previous case, we can guess that solutions of the zeroth order
equation are of the form $z^{i\alpha}$ and that $\alpha$ satisfies the
dispersion relation associated with the model. In the comoving frame, the
eikonal equation for the model provides us with the following equation:

$\left(k^{2}-\frac{\omega^{2}}{c^{2}}\right)\frac{1}{\chi\omega_{0}^{2}}(\omega_{0}^{2}-\gamma^{2}(\omega+vk)^{2})-\frac{1}{c^{2}}\gamma^{2}(\omega+vk)^{2}=0.$
(7.14)

It is a quartic equation whose roots $k_{i}(\omega)$, $i=1,2,3,4$ cannot be
managed in simple formulas unless some kind of approximation is provided. Our
ansatz is the following. We choose $\delta\to 0$ as an expansion parameter and
we put

$\displaystyle k$ $\displaystyle=:\frac{1}{\delta}u-\frac{\omega}{v},$ (7.15)
$\displaystyle\omega_{0}$ $\displaystyle=:\frac{1}{\delta}\bar{\omega}_{0}.$
(7.16)

We mean to indicate an expansion where $k$, unless it is zero, dominates over
$\omega/v$, i.e. $\omega$ is small relative to a nonzero $k$, and that also
$\omega_{0}$ is big compared to $\omega$ (this kind of approach leads also to
the Cauchy approximation, cf. e.g. [22], ch. 9). As a consequence, we find the
following re-writing of (7.14):

$\left(u^{2}-2\delta\frac{\omega}{v}u+\delta^{2}\frac{\omega^{2}}{\gamma^{2}v^{2}}\right)\left(1-\frac{\gamma^{2}v^{2}}{\bar{\omega}_{0}^{2}}u^{2}\right)-\chi\gamma^{2}v^{2}u^{2}=0.$
(7.17)

We look for a series solution in $u$:

$u=u_{0}+\delta\,u_{1}+\delta^{2}\,u_{2}+\delta^{3}\,u^{3}+\ldots$

As a consequence, we expect

$k=\frac{1}{\delta}\,u_{0}+\left(u_{1}-\frac{\omega}{v}\right)+\delta\,u_{2}+\delta^{2}\,u_{3}+\ldots$

We find at the zeroth order two vanishing degenerate solutions $u_{0}=0$, to
be associated with the modes $H,B$, and also

$\displaystyle u_{0P}$ $\displaystyle:=\frac{\bar{\omega}_{0}}{\gamma
v}\sqrt{1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}},$ (7.18) $\displaystyle u_{0N}$
$\displaystyle:=-u_{0P}.$ (7.19)

Corrections for the two non-degenerate solutions in the first order are

$\displaystyle u_{1P}$
$\displaystyle:=-\frac{\omega}{v}\chi\gamma^{2}\frac{v^{2}}{c^{2}}\frac{1}{1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}},$
(7.20) $\displaystyle u_{1N}$ $\displaystyle:=u_{1P}.$ (7.21)

Corrections at the first order for the two degenerate zero roots $u_{0}=0$
arise from the second-order contribution to the dispersion relation (as the
first order contribution vanishes identically for $u_{0}=0$): one finds

$\displaystyle u_{1H}$
$\displaystyle:=\frac{\omega}{v}\left(1+\sqrt{1+\chi}\frac{v}{c}\right)\frac{1}{1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}},$
(7.22) $\displaystyle u_{1B}$
$\displaystyle:=\frac{\omega}{v}\left(1-\sqrt{1+\chi}\frac{v}{c}\right)\frac{1}{1-\chi\gamma^{2}\frac{v^{2}}{c^{2}}}.$
(7.23)

We are not interested in further corrections. It is nice to point out that the
expressions we have found are compatible with the WKB behaviour of the
solutions found in [13].

###  7.3. Scattering coefficients

We proceed in the same way as exposed in Sections 4 and 5, performing a
perturbative expansion of Eq. (7.13) in the parameter $\eta$. We just give the
results of the scattering coefficients, analogous to (5.2):

$\displaystyle C_{1}$ $\displaystyle=\frac{1}{1+i\eta\tau}\,,$ (7.24)
$\displaystyle C_{2}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\frac{\pi(k_{P}^{2}-\omega^{2})(-1)^{i(k_{H}-k_{P})/(2\tilde{\beta})}}{2\tilde{\beta}^{2}V^{2}\epsilon^{2}(k_{P}-k_{N})(k_{P}-k_{B})}\text{csch}\left(\frac{\pi\gamma(k_{H}-k_{P})}{2\tilde{\beta}}\right)\,,$
(7.25) $\displaystyle C_{3}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\frac{\pi(k_{N}^{2}-\omega^{2})(-1)^{i(k_{H}-k_{N})/(2\tilde{\beta})}}{2\tilde{\beta}^{2}V^{2}\epsilon^{2}(k_{N}-k_{P})(k_{N}-k_{B})}\text{csch}\left(\frac{\pi\gamma(k_{H}-k_{N})}{2\tilde{\beta}}\right)\,,$
(7.26) $\displaystyle C_{4}$
$\displaystyle=\frac{i\eta}{1+i\eta\tau}\frac{\pi(k_{B}^{2}-\omega^{2})(-1)^{i(k_{H}-k_{B})/(2\tilde{\beta})}}{2\tilde{\beta}^{2}V^{2}\epsilon^{2}(k_{B}-k_{P})(k_{B}-k_{N})}\text{csch}\left(\frac{\pi\gamma(k_{H}-k_{B})}{2\tilde{\beta}}\right)\,,$
(7.27) $\displaystyle\tau$
$\displaystyle=\frac{(k_{H}^{2}-\omega^{2})}{\tilde{\beta}V^{2}\gamma\epsilon^{2}(k_{H}-k_{P})(k_{H}-k_{N})(k_{H}-k_{B})}\,.$
(7.28)

Recalling that, from (7.14),

$\displaystyle(k_{J}^{2}-\omega^{2})=\frac{\chi\omega_{0}^{2}\,\gamma^{2}\,(\omega+Vk_{J})^{2}}{\omega_{0}^{2}-\gamma^{2}(\omega+Vk_{J})^{2}}\overset{\omega_{lab}^{2}\ll\omega_{0}^{2}}{\approx}\chi\gamma^{2}(\omega+Vk_{J})^{2}\,,$
(7.29)

we can see that the expressions of the coefficients (7.24) basically reduce to
(5.2) at low frequencies. This fact was expected since the Cauchy dispersion
relation in an approximation of the Sellmeier for
$\omega_{lab}^{2}\ll\omega_{0}^{2}$, but it can be viewed also as a check of
consistency and robustness of our results.

##  8\. Comparison with the Orr-Sommerfeld approach

Equations like (3.8) are called of generalized Orr-Sommerfeld type. Such
equations emerge often in Analogue Gravity and have been studied extensively
in [13]. Here the authors developed a general technique for computing the
Hawking spectrum in the _transcritical_ case, using a perturbative approach in
the low-dispersion parameter. In this Section we compare our results to those
derived with the Orr-Sommerfeld approach: in particular, we establish a
relation between the effective temperature we have defined in Eq. (6.15) for
the subcritical regime, with the Hawking temperature the authors find in [13]
for the transcritical case. From equations (2.5) and (3.4) we we can identify
the low-dispersion limit as

$\displaystyle\epsilon^{2}:=\frac{1}{g^{2}}\rightarrow
0\,,\quad\frac{\mu^{2}}{g^{2}}=:\theta=\text{const}\,,\quad\frac{\eta}{g^{2}}=:\rho=\text{const}\,.$
(8.1)

Notice that this means that $\mu^{2}\propto\epsilon^{-2}$ and
$\eta\propto\epsilon^{-2}$. By eliminating the third order term, Eq. (3.8) can
be written as

$\displaystyle\epsilon^{2}f^{(4)}(x)+p_{3}(x)f^{\prime\prime}(x)+p_{2}(x)f^{\prime}(x)+p_{1}(x)f(x)=0\,,$
(8.2)

where

$\displaystyle p_{3}(x)$ $\displaystyle=$
$\displaystyle\frac{\left(g^{2}-\mu^{2}V^{2}+\eta
V^{2}\text{sech}^{2}(\tilde{\beta}x)\right)}{g^{2}V^{4}\gamma^{2}}=\frac{\left(1-\theta
V^{2}+\rho V^{2}\text{sech}^{2}(\tilde{\beta}x)\right)}{V^{4}\gamma^{2}}$
(8.3) $\displaystyle p_{2}(x)$ $\displaystyle=$
$\displaystyle\frac{2i\left(g^{2}\left(V^{2}-1\right)\omega-4\tilde{\beta}\eta
V^{3}\tanh(\tilde{\beta}x)\text{sech}^{2}(\tilde{\beta}x)\right)}{g^{2}V^{5}\gamma^{2}}$
(8.4) $\displaystyle=$
$\displaystyle\frac{2i\left(\left(V^{2}-1\right)\omega-4\tilde{\beta}\rho
V^{3}\tanh(\tilde{\beta}x)\text{sech}^{2}(\tilde{\beta}x)\right)}{V^{5}\gamma^{2}}\,,$
(8.5) $\displaystyle p_{1}(x)$ $\displaystyle=$
$\displaystyle\frac{\left(-2\tilde{\beta}^{2}\eta
V^{4}\text{sech}^{4}(\tilde{\beta}x)+4\tilde{\beta}^{2}\eta
V^{4}\tanh^{2}(\tilde{\beta}x)\text{sech}^{2}(\tilde{\beta}x)-g^{2}\left(V^{2}-1\right)^{2}\omega^{2}\right)}{g^{2}V^{6}\gamma^{2}}$
(8.6) $\displaystyle=$ $\displaystyle\frac{\left(-2\tilde{\beta}^{2}\rho
V^{4}\text{sech}^{4}(\tilde{\beta}x)+4\tilde{\beta}^{2}\rho
V^{4}\tanh^{2}(\tilde{\beta}x)\text{sech}^{2}(\tilde{\beta}x)-\left(V^{2}-1\right)^{2}\omega^{2}\right)}{V^{6}\gamma^{2}}\,.$
(8.7)

Notice that the functions $p_{j}(x)$ do not depend on $\epsilon$, as in the
Corley model. This form is directly comparable with the form in [13]. The
position of the horizon is defined by

$\displaystyle p_{3}(x)=0\,.$ (8.8)

By changing variable to $z=-e^{2\tilde{\beta}x}$ as before, we find

$\displaystyle z_{H\pm}$ $\displaystyle=\frac{g^{2}+2\eta
V^{2}-\mu^{2}V^{2}\pm 2\sqrt{g^{2}\eta
V^{2}+\eta^{2}V^{4}-\eta\mu^{2}V^{4}}}{g^{2}-\mu^{2}V^{2}}$ (8.9)
$\displaystyle=\frac{1+2\rho V^{2}-\theta V^{2}\pm 2\sqrt{\rho
V^{2}+\rho^{2}V^{4}-\rho\theta V^{4}}}{1-\theta V^{2}}\,.$ (8.10)

The corresponding points $x_{H\pm}$ are real if

$\displaystyle\Big{\\{}\eta<-\frac{(g^{2}-\mu^{2}V^{2})}{V^{2}}=:-\eta_{\text{min}}\Big{\\}}\cup\Big{\\{}\eta>0\Big{\\}}\,.$
(8.11)

In this comparison with the critical regime, we just consider the case of
positive Kerr nonlinearity: $\delta n>0\implies\eta<0$. Thus we neglect
positive values of $\eta$, and we identify $\eta<-\eta_{\text{min}}$ as the
condition of transcriticality of the perturbation. Notice that
$\eta_{\text{min}}>0$. In Figure 6 we plot the dispersion relation in the
critical case.

Figure 6.  The dispersion relation (3.5) represented at the peak of the
background field, for $\eta<-\eta_{\text{min}}$. In this configuration, for
any $\omega$, the modes $H$ and $P$ become imaginary as they experience an
event horizon: this is referred to as _transcritical regime_.

In [13] thermal Hawking radiation was predicted with a temperature

$\displaystyle
T_{H}=\frac{\kappa}{2\pi}=\frac{\gamma^{2}V^{2}n^{\prime}(x_{H})}{2\pi}\,,$
(8.12)

where $n^{\prime}(x_{H})$ is the derivative of the refraction index at the
horizon. We compute $n^{\prime}$ from the Cauchy dispersion formula including
the background correction in the comoving frame:

$n^{2}(x)=n_{0}^{2}-\frac{\lambda}{2g^{2}}\psi^{2}_{B}(x),$ (8.13)

where

$n_{0}^{2}:=\frac{\mu^{2}}{g^{2}}\,.$ (8.14)

Here we do not consider the dependence on $\omega_{lab}$ of the refractive
index (Eq. (3.5)), because the Orr-Sommerfeld approach is true in the low-
dispersion limit ($\omega_{lab}\sim 0$): indeed, the event horizon is defined
by

$\displaystyle p_{3}(x_{H})=0\quad\iff\quad n^{2}(x_{H})=\frac{1}{V^{2}}\,.$
(8.15)

With these considerations, we find

$\displaystyle
n^{\prime}(x_{H})=\frac{V}{2}\left(n^{2}(x_{H})\right)^{\prime}=V\left(-\frac{\lambda}{2g^{2}}\psi_{B}^{2}(x_{H})\right)\frac{\psi_{B}^{\prime}(x_{H})}{\psi_{B}(x_{H})}\,,$
(8.16)

and using (8.13), (8.15) and (3.1), we and up with

$\displaystyle T_{H}=$ $\displaystyle\frac{\beta\gamma
V}{2\pi}\frac{g^{2}-\mu^{2}V^{2}}{g^{2}}\frac{1+z_{H}}{1-z_{H}}$ (8.17)
$\displaystyle=$ $\displaystyle\frac{\beta\gamma
V}{2\pi}\frac{g^{2}-\mu^{2}V^{2}}{g^{2}}\sqrt{\frac{|\eta|-\eta_{\text{min}}}{|\eta|}}\,.$
(8.18)

This result is surprisingly similar to (6.15), although they were computed in
very different ways, one within the critical regime and the other in the
subcritical. We notice that if the expression (6.15) for $T_{pert}$ is
computed for near critical velocity

$n_{0}^{2}=\frac{\mu^{2}}{g^{2}}\lesssim\frac{1}{V^{2}}\,,$

and compare it to (8.17) in the very critical case
($|\eta|\gg\eta_{\text{min}}$), we find

$\displaystyle T_{pert}\approx\frac{T_{H}}{3}\,.$ (8.19)

We are not able to give an interpretation of the missing factor $3$ with
respect to the case of real turning points. We may suggest that it could be an
effect of considering a non-monotonic background such as (2.11), whereas the
results of the Orr-Sommerfeld approach in [13] were computed assuming a
monotonic background. In any case, the similarity of the critical and
subcritical temperatures is interesting and it suggests that they are related
by some physical mechanism: this will be a matter for future studies.

The same estimations can be done using the $\phi\psi$ model, introduced in
Section 7. We set $c=1$. The temperature estimated from the coefficients
(7.24) is

$\displaystyle\beta_{pert}=\frac{1}{T_{pert}}=\frac{2\pi\gamma}{2\beta}\lim\limits_{\omega\rightarrow
0}\frac{2k_{H}-k_{P}-k_{N}}{\omega}=\frac{2\pi\left(1+\gamma^{2}V^{2}\chi+V\sqrt{\chi+1}\right)}{\beta\gamma
V\left(1-V^{2}(\chi+1)\right)}\,.$ (8.20)

For near-critical velocity

$\displaystyle n_{0}^{2}=1+\chi\lesssim\frac{1}{V^{2}}\,,$ (8.21)

we see that the term between brackets in the numerator of (8.20) is

$\displaystyle\left(1+\gamma^{2}V^{2}\chi+V\sqrt{\chi+1}\right)\approx 3\,,$
(8.22)

which gives

$\displaystyle T_{pert}\approx\frac{\beta\gamma
V\left(1-V^{2}(\chi+1)\right)}{6\pi}\,.$ (8.23)

The analysis of the critical case with the Orr-Sommerfeld approach gives that
criticality is found for

$\displaystyle\eta<-\eta_{\text{min}}\,,\quad\eta_{\text{min}}:=\frac{1}{\chi}(1-\chi\gamma^{2}V^{2})>0\,$

and the Hawking temperature is

$\displaystyle T_{H}=$ $\displaystyle\frac{\beta}{2\pi\chi
V\gamma}(1-V^{2}(1+\chi))\sqrt{\frac{|\eta|-\eta_{\text{min}}}{|\eta|}}$
(8.24) $\displaystyle\approx$ $\displaystyle\frac{\beta\gamma
V}{2\pi}(1-V^{2}(1+\chi))\sqrt{\frac{|\eta|-\eta_{\text{min}}}{|\eta|}}\,,$
(8.25)

where in the last line we used (8.21) as before. We see that, once again, the
relation between the two temperatures (for $|\eta|\gg\eta_{\text{min}}$) is

$\displaystyle T_{pert}\approx\frac{T_{H}}{3}\,,$ (8.26)

confirming what we found using the Cauchy model.

##  9\. Conclusions

We proposed a new approach to Hawking-like radiation in the subcritical case
for a particular, but arguably realistic, class of soliton-like backgrounds in
a nonlinear dielectric. The method allows a straightforward analytical
solution of the scattering problem at the leading perturbation order of the
amplitude of the background field, represented by the parameter $\eta$. With
respect to other existing techniques, our approach does not rely on the
approximation of weak dispersion, and indeed our predictions are not
restricted to the $\omega\sim 0$ region. We tested our approach on a
simplified model of scalar electrodynamics and checked its robustness on the
scalar field reduction of the Hopfield model, which strongly corroborates the
results obtained in the simplified model. In both cases, we can define an
effective temperature associated with the spectrum of the emitted radiation,
which, in the limit where the subcritical case approaches the transcritical
one, is one-third of the Hawking temperature estimated by other established
methods for the transcritical regime. The interpretation of this fact is not
yet available to us: future work will be focused on the transition to the
critical regime, and it will possibly give us a clearer view.

The perturbative expansion we propose, although being naturally suited for the
study of the subcritical regime, is not theoretically limited to this case:
the problem of the transition to the transcritical regime is configured, in
this context, as possibly a matter of being able to compute enough
perturbative orders. This can be corroborated by the fact that a new
mathematical perspective on the phenomenon of particle creation was provided,
as we related particle creation to the solution of a fourth-order Fuchsian
equation. Fuchsian equations have many well-studied properties and much of
that theoretical machinery may be applied to the problem of analogue Hawking
radiation, potentially extending our analytical comprehension of the
phenomenon. This will be a matter for future studies.

##  Appendix A The conserved current

We compute the current vector, evaluated on normal modes $e^{ik_{\mu}x^{\mu}}$
in a generic frame moving with velocity $V$ with respect to the laboratory.
From Eq. (3.2) we find that if $\psi(x)=e^{ik_{\mu}x^{\mu}}$, then
$\phi(x)=-i\frac{gn^{\mu}k_{\mu}}{(v^{\mu}k_{\mu})^{2}}e^{ik_{\mu}x^{\mu}}$:
substituting these fields into (2.7) we find

$\displaystyle J^{\mu}\propto
v^{\mu}\,\omega_{lab}\left(1+\frac{\mu^{2}+\omega_{lab}^{2}}{\omega_{lab}^{2}}\right)+n^{\mu}\,\frac{g^{2}k_{lab}}{\omega_{lab}^{2}}\,,$
(A.1)

where $v^{\mu}=(\gamma,-V\gamma)$, $n^{\mu}=(-\gamma V,\gamma)$,
$\omega_{lab}=v^{\mu}k_{\mu}=\gamma(\omega+Vk)$,
$k_{lab}=-n^{\mu}k_{\mu}=\gamma(k+V\omega)$. We can verify that $J^{\mu}$ is
time-like. Indeed,

$\displaystyle
J^{\mu}J_{\mu}=\omega_{lab}^{2}\left(1+\frac{2(\mu^{2}+\omega_{lab}^{2})}{\omega_{lab}^{2}}+\frac{(\mu^{2}+\omega_{lab}^{2})^{2}}{\omega_{lab}^{4}}-\frac{g^{2}(\mu^{2}+\omega_{lab}^{2})}{\omega_{lab}^{4}}\right)\,,$
(A.2)

which is positive provided that $\mu^{2}>g^{2}$, which is the condition we
required at the beginning (see below Eq. (2.5)). Thus $J^{\mu}$ is time-like,
so the sign of $J^{0}$ is constant in every inertial reference frame: since in
the lab frame $\text{sign}(J^{0})=\text{sign}(\omega_{lab})$, this must be
true in every frame.

We now define the dispersion relation function

$\displaystyle\text{DR}(\omega,k):=\mu^{2}+\omega_{lab}^{2}-\frac{g^{2}k_{lab}^{2}}{\omega_{lab}^{2}}\,.$
(A.3)

The free normal modes of the theory satisfy $DR(\omega,k)=0$. It is easy to
show that the following relations hold

$\displaystyle J^{0}$
$\displaystyle\propto\partial_{\omega}\text{DR}|_{\text{DR}=0}\,,$ (A.4)
$\displaystyle J^{x}$
$\displaystyle\propto-\partial_{k}\text{DR}|_{\text{DR}=0}\,.$ (A.5)

Thus, we can identify the measure of the Hilbert space of free normal modes as

$\displaystyle
d\mu_{j}:=\frac{d\omega}{2\pi\partial_{k}\text{DR}|_{k=k_{j}(\omega)}}\,,$
(A.6)

where we are writing the field theory in the frequency representation and
$k_{j}$ are the different real solutions of $\text{DR}=0$. It is also
straightforward to show

$\displaystyle J^{x}=v_{g}(\omega)\,J^{0}\,,$ (A.7)

where $v_{g}(\omega)=\left.\frac{d\omega}{dk}\right|_{\text{DR}=0}$ is the
group velocity of the normal modes. This relation could be derived also from
the implicit function theorem, which states
$\frac{\partial_{k}\text{DR}}{\partial_{\omega}\text{DR}}=-\left.\frac{d\omega}{dk}\right|_{\text{DR}=0}$.

We can write the flux $J^{x}$ also in a more convenient way for future
computations. Since we called $k_{H}(\omega)$, $k_{P}(\omega)$,
$k_{N}(\omega)$ and $k_{B}(\omega)$ the four solutions of the dispersion
relation, we can write Eq. (A.3) as

$\displaystyle\text{DR}(\omega,k)=\frac{\gamma^{2}V^{4}}{\omega_{lab}^{2}}\prod_{j=H,P,N,B}(k-k_{j}(\omega))\,.$
(A.8)

It is now easy t verify that

$\displaystyle\partial_{\omega}\text{DR}|_{k=k_{i}(\omega)}=-\frac{\gamma^{2}V^{4}}{\omega_{lab}^{2}|_{k_{i}}}\frac{\partial
k_{i}}{\partial\omega}\prod_{j\neq i}(k_{i}(\omega)-k_{j}(\omega))\,,$ (A.9)

and then

$\displaystyle J^{x}(e^{-i\omega
t+ik_{i}(\omega)x})=v_{i}(\omega)\partial_{\omega}\text{DR}|_{k_{i}}=-\frac{\gamma^{2}V^{4}}{\omega_{lab}^{2}|_{k_{i}}}\prod_{j\neq
i}(k_{i}(\omega)-k_{j}(\omega))\,.$ (A.10)

##  Appendix B The conservation law

To understand the conservation law (6.5), we must consider that the initial
and final states of the scattering should be wave packets rather than plane
waves. A normalized wave packet centred around the frequency $\omega$ and
momentum $k_{j}(\omega)$ (here $k_{j}(\omega)$ is any of the solutions of the
dispersion relation) has the form

$\displaystyle\int\frac{d\omega}{2\pi\partial_{k}DR|_{k_{j}(\omega)}}f_{\epsilon}(\omega)\,e^{i\omega
t-ik_{j}(\omega)x}\,,$ (B.1)

where $\frac{d\omega}{2\pi\partial_{k}DR|_{k_{j}(\omega)}}$ is the measure of
the Hilbert space defined by plane waves (in the frequency representation) and
$f_{\epsilon}(\omega)$ is a distribution centred around $\omega$ with a small
bandwidth $\epsilon$, which satisfies

$\displaystyle\int\frac{d\omega}{2\pi\partial_{k}DR|_{k_{j}(\omega)}}|f_{\epsilon}(\omega)|^{2}=1\,.$

We can see Eq. (B.1) represents a wave packet by evaluating it at $x=vt$, for
some velocity $v$: for $t=\pm\infty$ we can say that the dominant contribution
comes from the stationary point, which is

$1-V\frac{\partial k_{j}}{\partial\omega}=0\,,$

which means that the function (B.1) is travelling at a constant velocity
$v=v_{j}(\omega)$, where $v_{j}(\omega)$ is the group velocity of the mode
$k_{j}$.

The asymptotic plane-wave solution (5.1) corresponds to the wave packets

$\displaystyle\psi(x,t)\sim\int\frac{d\omega}{2\pi\partial_{k}DR|_{k_{H}(\omega)}}f_{\epsilon}(\omega)\,e^{i\omega
t}\begin{cases}e^{-ik_{H}x}+C_{2}e^{-ik_{P}x}+C_{3}e^{-ik_{N}x}+C_{4}e^{-ik_{B}x}\,,&x\rightarrow-\infty\,,\\\
C_{1}e^{-ik_{H}x}\,,&x\rightarrow+\infty\,.\end{cases}$ (B.2)

At $t=-\infty$ the only contribution comes from the in-going $H$-mode, so the
measure of integration is chosen such that

$\left(\psi|\psi\right)_{t=-\infty}=1\equiv\int dxJ^{0}(\psi(x,t=-\infty))\,.$

The current conservation implies

$\displaystyle\frac{d}{dt}\int dxJ^{0}=\int dx\partial_{0}J^{0}=\int
dx\partial_{x}J^{x}=0\,.$ (B.3)

At $t=+\infty$ the solution splits into four outgoing localized wave packets
with central momenta: each mode propagates with a different group velocity, so
we can assume that the wave packets are not overlapping at $t=+\infty$. With
these considerations, the computation of the norm gives

$\displaystyle\int\frac{d\omega}{2\pi\partial_{k}\text{DR}|_{k_{H}(\omega)}}\frac{d\omega^{\prime}}{2\pi\partial_{k}\text{DR}|_{k_{H}(\omega^{\prime})}}f_{\epsilon}(\omega)f_{\epsilon}^{*}(\omega^{\prime})\sum_{i=H,P,N,B}|C_{i}|^{2}2\pi\partial_{\omega}\text{DR}|_{k_{i}}\delta(k_{H}(\omega)-k_{H}(\omega^{\prime}))=$
(B.4)
$\displaystyle=\int\frac{d\omega}{2\pi\partial_{k}\text{DR}_{k_{H}}}|f_{\epsilon}(\omega)|^{2}\frac{1}{2\pi\partial_{k}\text{DR}_{k_{H}}}\sum_{i=H,P,N,B}|C_{i}|^{2}2\pi
v_{i}\partial_{\omega}\text{DR}|_{k_{i}}\equiv 1\,.$ (B.5)

The last equality comes from current conservation. By recalling
$\partial_{k}DR|_{k_{j}(\omega)}=-v_{j}(\omega)\partial_{\omega}DR|_{k_{j}(\omega)}$
(from the implicit function theorem) we see that the sign of each term is
determined by the product $v_{j}(\omega)\,\omega_{lab}$. Since the equality
must hold for every distribution $f_{\epsilon}$, we finally obtain the
following relation between the absolute values,

$\displaystyle
1=|C_{1}|^{2}+|C_{2}|^{2}\left|\frac{v_{P}\partial_{\omega}DR|_{k_{P}}}{v_{H}\partial_{\omega}DR|_{k_{H}}}\right|-|C_{3}|^{2}\left|\frac{v_{N}\partial_{\omega}DR|_{k_{N}}}{v_{H}\partial_{\omega}DR|_{k_{H}}}\right|+|C_{4}|^{2}\left|\frac{v_{B}\partial_{\omega}DR|_{k_{B}}}{v_{H}\partial_{\omega}DR|_{k_{H}}}\right|\,,$
(B.6)

which is precisely (6.5).

##  Appendix C Solution of first order equation

We start computing the first-order solution $u_{1}(z)$ (Eq. (4.9)) around
$z=0$ ($x=-\infty$). We proceed as exposed in Section 4 to derive the
expression of the first-order solution. In the scattering of the $H$-mode, we
expect the asymptotic coefficient of the $H$-mode to be $O(1)$, while all the
others should be $O(\eta)$ (see Section 5): for this reason, the zero-order
solution we put into the source term (4.8) corresponds to a single $H$-mode,
which means

$\displaystyle u_{0}(z)=z^{i\alpha_{H}}\,,\quad
i\alpha_{H}=i\frac{k_{H}}{2\beta}+\frac{3}{2}+i\frac{\Omega}{V}\,.$ (C.1)

From this choice, we get

$\displaystyle r_{1}(z)$
$\displaystyle=\frac{z^{-3+i\alpha_{H}}\left(4\alpha_{H}^{2}(z-1)^{2}+4i\alpha_{H}\left(5z^{2}-6z+1\right)-25z^{2}+2z-1\right)}{4(z-1)^{4}V^{2}y^{2}}\,,$
(C.2) $\displaystyle\frac{W_{1}(z)}{W(z)}$
$\displaystyle=\frac{i\left(4\alpha_{H}^{2}(z-1)^{2}+4i\alpha_{H}\left(5z^{2}-6z+1\right)-25z^{2}+2z-1\right)}{4(z-1)^{4}V^{2}\gamma^{2}(\alpha_{H}-\alpha_{P})(\alpha_{H}-\alpha_{N})(\alpha_{H}-\alpha_{B})}\,,$
(C.3) $\displaystyle\frac{W_{2}(z)}{W(z)}$
$\displaystyle=-\frac{i\left(4\alpha_{H}^{2}(z-1)^{2}+4i\alpha_{H}\left(5z^{2}-6z+1\right)-25z^{2}+2z-1\right)z^{i(\alpha_{H}-\alpha_{P})}}{4(z-1)^{4}V^{2}\gamma^{2}(\alpha_{P}-\alpha_{H})(\alpha_{P}-\alpha_{N})(\alpha_{P}-\alpha_{B})}\,,$
(C.4) $\displaystyle\frac{W_{3}(z)}{W(z)}$
$\displaystyle=-\frac{i\left(4\alpha_{H}^{2}(z-1)^{2}+4i\alpha_{H}\left(5z^{2}-6z+1\right)-25z^{2}+2z-1\right)z^{i(\alpha_{H}-\alpha_{N})}}{4(z-1)^{4}V^{2}\gamma^{2}(\alpha_{N}-\alpha_{H})(\alpha_{N}-\alpha_{P})(\alpha_{N}-\alpha_{B})}\,,$
(C.5) $\displaystyle\frac{W_{4}(z)}{W(z)}$
$\displaystyle=-\frac{i\left(4\alpha_{H}^{2}(z-1)^{2}+4i\alpha_{H}\left(5z^{2}-6z+1\right)-25z^{2}+2z-1\right)z^{i(\alpha_{H}-\alpha_{B})}}{4(z-1)^{4}V^{2}\gamma^{2}(\alpha_{B}-\alpha_{H})(\alpha_{B}-\alpha_{P})(\alpha_{B}-\alpha_{N})}\,.$
(C.6)

Computing the integrals in (4.9), we get

$\displaystyle\tilde{u}_{1}(z)$
$\displaystyle=\frac{i\,z^{i\alpha_{H}}}{4V^{2}\gamma^{2}}\left(\frac{-4\alpha_{H}^{2}(z-1)^{2}-4i\alpha_{H}\left(5z^{2}-8z+3\right)+25z^{2}-26z+9}{(z-1)^{3}(\alpha_{H}-\alpha_{P})(\alpha_{H}-\alpha_{N})(\alpha_{H}-\alpha_{B})}+\right.$
(C.7)
$\displaystyle\left.\frac{z\left(i(2\alpha_{H}+5i)^{2}\,_{2}F_{1}(2,i(\alpha_{H}-\alpha_{P})+1,i(\alpha_{H}-\alpha_{P})+2;z)+16(\alpha_{H}+3i)\,_{2}F_{1}(3,i(\alpha_{H}-\alpha_{P})+1,i(\alpha_{H}-\alpha_{P})+2;z)\right)}{(\alpha_{H}-\alpha_{P})(\alpha_{H}-\alpha_{P}-i)(\alpha_{P}-\alpha_{N})(\alpha_{P}-\alpha_{B})}\right.$
(C.8)
$\displaystyle\left.-\frac{z\left(24\,i\,_{2}F_{1}(4,i(\alpha_{H}-\alpha_{P})+1,i(\alpha_{H}-\alpha_{P})+2;z)\right)}{(\alpha_{H}-\alpha_{P})(\alpha_{H}-\alpha_{P}-i)(\alpha_{P}-\alpha_{N})(\alpha_{P}-\alpha_{B})}+\right.$
(C.9)
$\displaystyle\left.\frac{z\left(i(2\alpha_{H}+5i)^{2}\,_{2}F_{1}(2,i(\alpha_{H}-\alpha_{N})+1,i(\alpha_{H}-\alpha_{N})+2;z)+16(\alpha_{H}+3i)\,_{2}F_{1}(3,i(\alpha_{H}-\alpha_{N})+1,i(\alpha_{H}-\alpha_{N})+2;z)\right)}{(\alpha_{H}-\alpha_{N})(\alpha_{H}-\alpha_{N}-i)(\alpha_{N}-\alpha_{P})(\alpha_{N}-\alpha_{B})}\right.$
(C.10)
$\displaystyle\left.-\frac{z\left(24\,i\,_{2}F_{1}(4,i(\alpha_{H}-\alpha_{N})+1,i(\alpha_{H}-\alpha_{N})+2;z)\right)}{(\alpha_{H}-\alpha_{N})(\alpha_{H}-\alpha_{N}-i)(\alpha_{N}-\alpha_{P})(\alpha_{N}-\alpha_{B})}+\right.$
(C.11)
$\displaystyle\left.\frac{z\left(i(2\alpha_{H}+5i)^{2}\,_{2}F_{1}(2,i(\alpha_{H}-\alpha_{B})+1,i(\alpha_{H}-\alpha_{B})+2;z)+16(\alpha_{H}+3i)\,_{2}F_{1}(3,i(\alpha_{H}-\alpha_{B})+1,i(\alpha_{H}-\alpha_{B})+2;z)\right)}{(\alpha_{H}-\alpha_{B})(\alpha_{H}-\alpha_{B}-i)(\alpha_{B}-\alpha_{P})(\alpha_{B}-\alpha_{N})}\right.$
(C.12)
$\displaystyle\left.-\frac{z\left(24\,i\,_{2}F_{1}(4,i(\alpha_{H}-\alpha_{B})+1,i(\alpha_{H}-\alpha_{B})+2;z)\right)}{(\alpha_{H}-\alpha_{B})(\alpha_{H}-\alpha_{B}-i)(\alpha_{B}-\alpha_{P})(\alpha_{B}-\alpha_{N})}\right)\,.$
(C.13)

This is a _particular_ solution of (4.5): the general expression is obtained
by adding a linear combination of free-field solutions (with coefficients of
order $O(\eta)$). The coefficients of such combinations will be determined
later based on the boundary conditions. From (C.7) we can compute the
asymptotic expression of $\tilde{u}_{1}$ around $z=0$ ($x=-\infty$):

$\displaystyle\tilde{u}_{1}(z\approx 0)$
$\displaystyle=c_{H}\,z^{i\alpha_{H}}\left(1+O(z)\right)\,,$ (C.14)
$\displaystyle c_{H}$
$\displaystyle=-\frac{i\left(-4\alpha_{H}^{2}-12i\alpha_{H}+9\right)}{4V^{2}\gamma^{2}(\alpha_{H}-\alpha_{P})(\alpha_{H}-\alpha_{N})(\alpha_{H}-\alpha_{B})}\,.$
(C.15)

The function (C.7) is defined on the whole complex plane: to write the
asymptotic expression at $z=\infty$ ($x=+\infty$) we use the connection
formulas of the hypergeometric function:

${}_{2}F_{1}(a,b,c;z)=$
$\displaystyle\frac{\left(-z\right)^{-a}\Gamma(c)\Gamma(b-a)}{\Gamma(b)\Gamma(c-a)}\,_{2}F_{1}\left(a,a-c+1;a-b+1;\frac{1}{z}\right)$
(C.16)
$\displaystyle+\frac{\left(-z\right)^{-b}\Gamma(c)\Gamma(a-b)}{\Gamma(a)\Gamma(c-b)}\,_{2}F_{1}\left(b,b-c+1;-a+b+1;\frac{1}{z}\right)\,.$
(C.17)

Expanding around $z=\infty$ we find

$\displaystyle\tilde{u}_{1}(z\approx\infty)=$ $\displaystyle-
c_{P}\,z^{i\alpha_{P}}\left(1+O(z)\right)-c_{N}\,z^{i\alpha_{N}}\left(1+O(z)\right)-c_{B}\,z^{i\alpha_{B}}\left(1+O(z)\right)\,,$
(C.18) $\displaystyle c_{P}$
$\displaystyle=-\frac{i\pi(2\alpha_{P}+3i)^{2}\left(-1\right)^{i(\alpha_{H}-\alpha_{P})}\text{csch}(\pi(\alpha_{H}-\alpha_{P}))}{4V^{2}\gamma^{2}(\alpha_{P}-\alpha_{N})(\alpha_{P}-\alpha_{B})}\,,$
(C.19) $\displaystyle c_{N}$
$\displaystyle=-\frac{i\pi(2\alpha_{N}+3i)^{2}\left(-1\right)^{i(\alpha_{H}-\alpha_{N})}\text{csch}(\pi(\alpha_{H}-\alpha_{N}))}{4V^{2}\gamma^{2}(\alpha_{N}-\alpha_{P})(\alpha_{N}-\alpha_{B})}\,,$
(C.20) $\displaystyle c_{B}$
$\displaystyle=-\frac{i\pi(2\alpha_{B}+3i)^{2}\left(-1\right)^{i(\alpha_{H}-\alpha_{B})}\text{csch}(\pi(\alpha_{H}-\alpha_{B}))}{4V^{2}\gamma^{2}(\alpha_{B}-\alpha_{N})(\alpha_{B}-\alpha_{P})}\,.$
(C.21)

From the boundary conditions (5.1) we see that we do not have the modes $P$,
$N$ and $B$ at right infinity; so we add to the particular solution the linear
combination

$c_{P}z^{i\alpha_{P}}+c_{N}z^{i\alpha_{N}}+c_{B}z^{i\alpha_{B}}\,.$

The first-order correction to the solution is thus

$\displaystyle
u_{1}(z)=\tilde{u}_{1}(z)+c_{P}z^{i\alpha_{P}}+c_{N}z^{i\alpha_{N}}+c_{B}z^{i\alpha_{B}}\,,$

such that $u_{1}(z\approx\infty)=0$. Finally, we can write down the solution
$u(z)$ at first perturbative order and its asymptotic behaviour:

$\displaystyle
u(z)\sim\begin{cases}(1+\eta\,c_{H})\,z^{i\alpha_{H}}+\eta\,c_{P}\,z^{i\alpha_{P}}+\eta\,c_{N}\,z^{i\alpha_{N}}+\eta\,c_{B}\,z^{i\alpha_{B}}\,,&z\rightarrow
0\,,\\\ z^{i\alpha_{H}}\,,&z\rightarrow\infty\,.\end{cases}$ (C.22)

The coefficients (5.2) are obtained dividing (C.22) by $(1+\eta\,c_{H})$ and
recalling $i\alpha_{j}=i\frac{k_{j}}{2\beta}+\frac{3}{2}+i\frac{\Omega}{V}$.

## References

  * [1] 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).
  * [2] S. Weinfurtner, E. W. Tedford, M. C. J. Penrice, W. G. Unruh and G. A. Lawrence, “Classical aspects of Hawking radiation verified in analogue gravity experiment,” Lect. Notes Phys. 870, 167 (2013).
  * [3] F. Michel and R. Parentani, “Probing the thermal character of analogue Hawking radiation for shallow water waves?,” Phys. Rev. D 90, no.4, 044033 (2014).
  * [4] L. P. Euvé, F. Michel, R. Parentani and G. Rousseaux, “Wave blocking and partial transmission in subcritical flows over an obstacle,” Phys. Rev. D 91, no. 2, 024020 (2015).
  * [5] S. Robertson, F. Michel and R. Parentani, “Scattering of gravity waves in subcritical flows over an obstacle,” Phys. Rev. D 93, no. 12, 124060 (2016).
  * [6] A. Coutant and S. Weinfurtner, “The imprint of the analogue Hawking effect in subcritical flows,” Phys. Rev. D 94, no. 6, 064026 (2016).
  * [7] H. Bremmer, “The W.K.B. approximation as the first term of a geometric-optical series,” Comm. Pure Appl. Math. 4 105 (1951).
  * [8] H. Bremmer, “The propagation of electromagnetic waves through a stratified medium and its W.K.B. approximation for oblique incidence,” Physica 15, 593 (1949).
  * [9] A. Coutant and S. Weinfurtner, “Low-frequency analogue Hawking radiation: The Korteweg-de Vries model,” Phys. Rev. D 97, no. 2, 025005 (2018).
  * [10] A. Coutant and S. Weinfurtner, “Low frequency analogue Hawking radiation: The Bogoliubov-de Gennes model,” Phys. Rev. D 97, no. 2, 025006 (2018).
  * [11] S. Finazzi and I. Carusotto, “Quantum vacuum emission in a nonlinear optical medium illuminated by a strong laser pulse,” Phys. Rev. A 87, no.2, 023803 (2013).
  * [12] S. Finazzi and I. Carusotto, “Spontaneous quantum emission from analog white holes in a nonlinear optical medium,” Phys. Rev. A 89, no. 5, 053807 (2014).
  * [13] F. Belgiorno, S. L. Cacciatori and A. Viganò, “Analog Hawking effect: A master equation,” Phys. Rev. D 102, no. 10, 105003 (2020).
  * [14] F. Belgiorno, S. L. Cacciatori, A. Farahat and A. Viganò, “Analog Hawking effect: BEC and surface waves,” Phys. Rev. D 102, no. 10, 105004 (2020).
  * [15] T. Nishimoto, “On the Orr-Sommerfeld type equations, II Connection formulas,” Kŏdai Math. Sem. Rep. 29, (1978), 233.
  * [16] T. Nishimoto, “On the Orr-Sommerfeld type equations, I W.K.B. approximation,” Kŏdai Math. Sem. Rep. 24, (1972), 281.
  * [17] F. Belgiorno, S. L. Cacciatori and F. Dalla Piazza, “Hawking effect in dielectric media and the Hopfield model,” Phys. Rev. D 91, no.12, 124063 (2015).
  * [18] F. Belgiorno and S. L. Cacciatori, “Analogous Hawking Effect in Dielectric Media and Solitonic Solutions,” Universe 6, no.8, 127 (2020).
  * [19] J. J. Hopfield, “Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals,” Phys. Rev. 112 (1958), 1555–1567.
  * [20] A.V. Chueshev, V.V. Chueshev, “Variational Formulas of the Monodromy Group for a Third-Order Equation on a Compact Riemann Surface,” Journal of Siberian Federal University. Mathematics & Physics 2022, 15(3), 308–318.
  * [21] F. Belgiorno, S. L. Cacciatori, F. Dalla Piazza and M. Doronzo, “Hopfield-Kerr model and analogue black hole radiation in dielectrics,” Phys. Rev. D 96, no.9, 096024 (2017).
  * [22] F. Belgiorno, S. L. Cacciatori and D. Faccio, Hawking Radiation, World Scientific Publishing Company, Singapore (2018).