Numerical analysis of backreaction in acoustic black holes

0.25 MB·English

by  Roberto Balbinot

Numerical analysis of backreaction in acoustic black holes
Roberto Balbinot1, Alessandro Fabbri2, Serena Fagnocchi3,1 and Alessandro Nagar4
1 Dipartimento di Fisica, Universit`a di Bologna and INFN sez. di Bologna, Via Irnerio 46, 40126 Bologna Italy
2 Departamento de Fisica Teorica and IFIC, Centro Mixto Universidad de Valencia–CSIC.
Facultad de Fis`ica, Universidad de Valencia, 46100 Burjassot (Valencia), Spain
3 Centro Enrico Fermi, Compendio Viminale, 00184 Roma, Italy
4 Dipartimento di Fisica, Politecnico di Torino and INFN sez. di Torino,
Corso Duca degli Abruzzi 24, 10129 Torino, Italy.
(Dated: February 7, 2008)
Using methods of Quantum Field Theory in curved spacetime, the first order in ~ quantum
correctionstothemotionofafluidinanacousticblackholeconfigurationarenumericallycomputed.
These corrections arise from the non linear backreaction of the emitted phonons. Time dependent
(isolated system) and equilibrium configurations (hole in a sonic cavity) are both analyzed.
6
PACSnumbers: 04.62.+v,04.70.Dy,47.40.Ki
0
0
2
I. INTRODUCTION II. THE ACOUSTIC BLACK HOLE

n
a
An acoustic black hole is a region of a fluid where its
J
Black hole radiation predicted by Hawking in 1974 [1] motion is supersonic. Here sound can not escape up-
0
is one ofthe most spectacularresults of moderntheoret- streambeing draggedbythe fluid. The boundaryofthis
2
ical physics. region is formed by sonic points where the speed of the

1 fluidequalsthelocalspeedofsound. Thisistheacoustic
Even more surprising is the fact that this effect is not
v horizon. A simple device to establish a transonic flow is
3 peculiar of gravitational physics, but is also expected in aconvergingdivergingdeLavalnozzle[3,4]. Forstation-
8 many completely different contexts of condensed matter
ary free fluid flow the acoustic horizon occurs exactly at
0 physics[2,4,5]. Afluidundergoingsupersonicmotionis
the waist of the nozzle.
1 thesimplestexampleofwhatonecallsan“acousticblack
The basic equations describing the system at the clas-
0
hole”. For this configuration Unruh [2], using Hawking
6 sicallevelarethecontinuityandtheBernoulliequations.
arguments, predicted an emission of thermal phonons.
0 We assume a one-dimensional stationary flow, therefore
This emission affects the behaviour of the underlying
/ all the relevant quantities depend on z only, the spatial
c fluid because of the non linearity of the hydrodynami-
q coordinaterunningalongthe axisofthe de Lavalnozzle.
cal equations governing its motion.
- The continuity equation then reads
r
g Using methods borrowedfrom QuantumField Theory
A(z)ρ(z)v(z)=const=D , (1)
:
v incurvedspacetime,thisquantumbackreactionhasbeen
i studiedforthefirsttimein[6],wherethethefirstorderin whereAistheareaofthetransversesectionofthenozzle,
X
~ corrections to the classical hydrodynamical equations ρthefluiddensityandv thefluidvelocity. TheBernoulli
ar were given. Because of intrinsic mathematical difficul- equation, under the above hypothesis, gives
ties, the analysis was restricted to the region very close
tothe “sonichorizon”oftheacousticblackhole;i.e.,the v2
+µ(ρ)=0 , (2)
region where the fluid motion changes from subsonic to 2
supersonic. There, analytical expressions for the quan-
tum corrections to the density and velocity of the mean whereµ(ρ)istheenthalpy. Wehavefurtherassumedthe
flow have been provided. fluid to be homentropic and irrotational. The speed of
sound c is defined as
However,tohaveadetaileddescriptionthroughoutthe
dµ
entiresystemonehastoproceedwithnumerics. Thiswill c2 =ρ . (3)
dρ
be the aim of our present paper, which is organized as
follows: in Sec. II we outline the classical fluid configu-
For constant c (the case we consider) integration of (3)
ration which describes an acoustic black hole; the quan-
gives
tumbackreactionequationsarediscussedinSec.III,with
emphasis on the choice of quantum state in which the ρ
µ(ρ)=c2ln , (4)
phonons field has to be quantized. In Secs. IV and V we ρ
0
give the numerical estimates for the quantum correction
to the mean flow in two different cases: isolated system which inserted in Bernoulli equation yields
and system in equilibrium in a sonic cavity respectively.
−v2
Section VI contains the final discussion. ρ=ρ0e 2c2 , (5)
2
x 104
x 10−8 0
3.5
−1
3 m/s] −2
c
v [ −3
2.5 −4
−5
2m] 2 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
c
A [ 0.06
1.5
0.05
3]0.04
1 m
g/c0.03
ρ [0.02
0.5
0.01
0
0 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
0.02 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 z [cm]
z [cm]
FIG.2: Velocity(top)anddensity(bottom)fromEq.(6). The
vertical dashed line correspond to the location of the sonic
FIG.1: DependenceofthecrosssectionofthedeLavalnozzle
on the position z for a velocity field given by Eq. (6) and horizon (zH = 0.025 cm). The sonic black hole corresponds
depicted in Fig. 2. The vertical dashed line corresponds to to z<zH.
thelocation of thesonic horizon zH =0.025 cm.
v. Quantization of these modes leads to the conclusion
with ρ a constant. The assumedconstancyof the speed that in presence of a sonic horizon a thermal emissionof
0
ofsoundalsogivesthepressurepasp=c2ρ. Thevelocity phononsisexpected,incompleteanalogyofwhatHawk-
profile describing the acoustic black hole is chosen as ing found for gravitational black holes. The emission
temperature of the phonons is T =~k/(2πcκ ), where
H B
2 k is the surface gravity of the sonic horizon, defined as
v =c arctan[β(z−z )]−1 , (6)
H
π
(cid:26) (cid:27) 1 d
k = (c2−v2) . (8)
where z = zH denotes the position of the waist of the 2dn
(cid:12)zH
nozzle (the sonic horizon). In the laboratory frame the (cid:12)
fluid is moving from right to left, so v <0 and the sonic κB istheBoltzmannconstantandn(cid:12)(cid:12) isthe normaltothe
horizon occurs where v = −c. The constant D entering horizon. For the specific acoustic black hole model we
the continuity equation is determined by requiring the consider in this paper, TH =1.1598×10−5 ◦K.
fluid to be sonic at the waist; i.e.,
D =−cA ρ e−1/2 =−AHpH , (7) III. THE BACKRECTION EQUATIONS
H 0
c
where A is the area at the horizon, p the pressure The phonons quantum emission previously discussed
H H
and ρ =p e1/2/c2. Given this, the profile of the nozzle modifies the underlying fluid flow accordingto the back-
0 H
can be computed from Eq. (1) and is depicted in Fig. 1, reaction equations derived in Ref. [6], to which we refer
where we have used A = 10−8 cm2, β = 600 cm−1, forfurtherdetails. Foraone-dimensionalflow,theyread
H
p =2×106 Pa and c =250 m/s. These latter two are
H 1
typical values for liquid Helium. The profiles of density Aρ +∂ (Aρ v )=∂ (hT(2)i+vhT(2)i) ,
B z B B z c tz zz
and velocity are shown in Fig. 2, where the significant (cid:20) (cid:21)
(9)
range of z is [0,0.05] cm; the horizon lies at z = 0.025
H
cm and its locationis indicated by a verticaldashed line v2 1
A ψ˙ + B +µ(ρ ) = hT(2)i . (10)
in the figures. In the region z > zH the motion of the B 2 B 2
fluid is subsonic; the acoustic black hole is the region (cid:20) (cid:21)
z <zH. Here ρB and vB are the quantum corrected density and
AsshownbyUnruh[2],soundwavespropagatinginan velocity fields and ψB is the velocity potential; i.e.,
inhomogeneous fluid are described as a massless scalar ∂zψB = vB; the overdot stands for time derivative. The
(2)
field propagating in an effective curved spacetime de- hT i which drive the backreactionare the quantum ex-
ab
scribed by an “acoustic metric” which depends on ρ and pectation values of the pseudo energy momentum tensor
3
6e-07 0.0002
Hartle-Hawking
Unruh
4e-07
0.0001
2e-07
0
0
F2 G2
1e-08
-2e-07 8e-09 -0.0001
6e-09
-4e-07 4e-09
-0.0002
2e-09
-6e-07 0
0.01 0.02 0.03 0.04 0.05
0 0.01 0.02 0.03 0.04 0.05 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
z [cm] z [cm]
FIG. 3: Quantum sources in the Unruh state (solid lines) and in the Hartle-Hawking state (dashed line). To appreciate the
difference in F2 taken in the Hartle-Hawking and Unruh states, we show it in the inset. The G2 is the same for both states.
The sonic horizon is located at zH =0.025 cm.
quadratic in the phonons field. To evaluate (ρB,vB) up and ∆±± are functions which depend on the choice of
to O(~) terms, the r.h.s. of the backreaction equations the quantum state of the phonons field. For the Unruh
(9) and (10) needs just to be evaluated on the classical state:
background (ρ,v) of Sec. II.
Thequantumstateofthefieldinwhichtheexpectation ∆ ≡∆U = 0 , (15)
++ ++
valueshavetobecomputeddependsonthephysicalsitu-
~k2
acatipoinngonpehownaonntssrtaoddiaetsicorniblee.adFsortoanaitsimolaetveadrihaotlieo,nthofetehse- ∆−− ≡∆U−− = 48πc4 . (16)
underlying medium, i.e. ρ (t,z) and v (t,z). The ap-
B B
propriatequantumstateinthiscaseistheanalogueofthe For the Hartle-Hawking state instead
Unruhstate[7]. Incasethesystemismaintainedinther-
mal equilibrium with the surroundings (that is, putting ~k2
a sonic cavity in the subsonic asymptotic right region), ∆±± ≡∆H±H± = . (17)
48πc4
thequantumstateistheanalogueoftheHartle-Hawking
state [8], the thermal equilibrium state at T = T . In
H
this case the system remains stationary, i.e. ρ (z) and From Eqs. (15)–(17) it follows that, in the asymp-
B
v (z). totic subsonic region z → +∞, ∆U±± describes a flux
B
Neglecting backscattering of the phonons, hT(2)i can ofphononsatatemperatureTH,whereas∆H±H± describes
ab atwo-dimensionalgasofphononsatthermalequilibrium
be approximated with the Polyakovstress tensor [9, 10].
at the temperature T . To find first order in ~ correc-
Introducing for the sake of simplicity null coordinates H
tions to the classical sonic black hole fluid configuration
(ρ(z),v(z)) described in Eqs. (5) and (6) we write
dz
x± =c t∓ (11)
c±v
(cid:18) Z (cid:19)
ψ = ψ(z)+ǫψ (t,z) , (18)
B 1
the Polyakovstress tensor reads: ρ = ρ(z)+ǫρ (t,z) , (19)
B 1
~
hT±(2±)i = −12πC1/2C,−±1±/2+∆±± (12) with vB = ∂zψB and ǫ is a dimensionless expansion pa-
~ rameter [11]:
hT(2)i = C−1(lnC) , (13)
6π ,+−
~
ǫ= . (20)
where |D|A
H
ρc2−v2
C = (14) For our system ǫ=1.317×10−14.
c c2
4
x 10−7
1 6
0.8 4
0.6 2
m/s] 0.4 3cm] 0 t=0
v [c1 0.2 ρ [g/1−2
t=0
0
−4
t=T
−0.2 end
−6 t=T
end
−0.4
0.016 0.018 0.02 0.022 0.024 0.026 0.028 0.03 0.032 −8
0.016 0.018 0.02 0.022 0.024 0.026 0.028 0.03 0.032
z [cm]
z [cm]
FIG. 4: Unruh state: time evolution of the backreaction equations for cκt ≪ 1. Snapshots of the quantum correction to the
velocity v1 (left panel) and to thedensity(right panel) for an evolution timetend =0.2tmax≈2.09×10−2 µs (cκt=0.2). The
vertical dashed line corresponds to the location of the classical sonic horizon. The delay between one snapshot and the other
(between t=0 and t=tend) is ∆t≈1.74×10−3 µs.
The backreactionequations linearized in ǫ then become IV. UNRUH STATE
ǫ Aρ˙ +∂ [A(ρ v+ρv )]
1 z 1 1
As said before, in Ref. [6] the backreaction equations
(cid:26) (cid:27)
hT(2)i hT(2)i wereanalyticallysolvedjustforz ≈zH toallowaTaylor
=c2∂z (c−++v)2 − (c+−−v)2 ≡ǫF2 , (21) expansion of the sources up to linear terms. In this sec-
" # tion we compute the numerical solutionall overthe noz-
c2 hT(2)i zle. We finite-difference the systemofEqs.(23) and(24)
ǫ A ψ˙1+vv1+ ρ1 = ≡ǫG2 . (22) and solve it numerically in the time domain as an ini-
ρ 2
(cid:20) (cid:18) (cid:19)(cid:21) tial value problem. The equations are discretized on an
evenlyspacedgrid(0,z )withz =0.05cm. Followinga
Using the backgroundequations (1) and (2), satisfied by c c
standardconventioninnumericalfluidmechanics[13],we
ρ and v, the continuity equation can be rewritten as
haveusedastaggeredgrid,i.e. bothz =0andz =z are
c
v2v′ ρv′ F thought to lie on cell interfaces while the hydrodynam-
′ ′ ′′ 2
ρ˙1+vρ1+ c2 ρ1− v ψ1+ρψ1 = A , (23) ics quantitiesaredefined oncellcenters. As a result,the
firstpointofourcomputationaldomainisz =∆z/2and
1
whereas the Bernoulli equation is the last is z =z −∆z/2. We notice that ∆z is cho-
imax c
senso thatthe horizonis locatedata cell interface. The
ψ˙ +vψ′ + c2ρ = G2 , (24) reasonforthisisthat,eventhoughthesquarebracketon
1 1 ρ 1 A the r.h.s. of Eq. (21) is analytically regular at z = zH
(where v =−c), the presence of the combination (c+v)
with a prime indicating derivative with respect to z. at denominator in Eq. (23) can give problems (i.e., di-
The profiles for the quantum sources F and G are vision by zero) due to the discretization procedure. The
2 2
depictedinFig.3forthe Unruhstate(solidline)andfor use of a staggered grid bypasses this difficulty since the
the Hartle-Hawking state (dashed line). The difference sonicpointturnsouttobe alwaysdisplacedwithrespect
betweenthestatesisreflectedonF only(andisshownin to the grid points.
2
the inset in the left panel of the figure), while G , being
2
relatedto the trace anomalywhich is state-independent, Initial conditions are chosen so that at t=0 the solu-
is unchanged. One can note the appearance of a maxi- tion is the classical one; i.e., ρ (t = 0) = ψ (t = 0) = 0.
1 1
mum and a minimum in the region z ∈ [0.02,0.03] cm. Then the backreaction is switched on. As in Ref. [6],
Outside this range,F and G rapidly dropto zero. The since the quantum sources are computed only for the
2 2
analysis of Ref. [6], being limited to the regionvery near static classical background, the validity of the solution
to z , could not catch this non trivial structure. is limited by the condition cκt ≪ 1, where we have in-
H
5
x 10−6
3.5 3
Hartle−Hawking Hartle−Hawking
Unruh Unruh
3
2.5
2.5
2
2
1.5
1.5
v [cm/s]10.51 3ρ [g/cm]1 1
0.5
0
0
−0.5
−0.5
−1
−1.5 −1
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
z [cm] z [cm]
FIG. 5: Profile of the quantum correction to the velocity v1 (left panel) and to ρ1 (right panel) due to backreaction in the
Hartle-Hawking state (black lines), compared with the profile of v1 and ρ1 at t = tend in the Unruh state (blue lines). The
vertical dashed lines correspond to theposition of thesonic horizon at t=0.
troduced the constant κ as covering the numerical domain) is sufficient to be in the
convergence regime (see Appendix A for discussion).
dv
κ=c−2k =c−1 (25) In Fig. 4 we have snapshots of the time evolution of
dz the profiles of v (left panel) and ρ (right panel). For
(cid:12)zH 1 1
(cid:12) this particular computation, we have considered a to-
(cid:12)
with dimension [length]−1. For the s(cid:12)onic black hole con- tal evolution time t = 0.2 t . The initial and fi-
end max
sideredhere,the shorttime conditiondeterminesamax- nal snapshots are depicted in red and blue respectively.
imum evolution time (cκt = 1) of t = 0.104 µs; so The time delay between one snapshot and the other is
max
it is possible to extractonly informations about how the ∆t ≃ 1.74×10−3µs. The quantum corrected velocity is
backreaction starts. obtained as v = ∂ ψ , the derivative being computed
1 z 1
Before discussing our numerical results, we briefly de- directly from the numerical data by means of a second
scribethenumericalalgorithmsimplemented,furtherde- order finite-difference approximation.
tails can be found in Appendix A. The numerical solution confirms the near horizon be-
We are dealing with a system of Partial Differen- haviour obtained in Ref. [6]: the fluid slows down close
tial Equations (PDEs), where the equation for ρ is of to the horizon (v > 0, remember that v < 0 because
1 1
convection-diffusiontype,due to the parabolictermpro- the fluid flows from right to left), causing the horizon to
portional to ψ′′, while the equation for ψ is a simple movetotheleft,andthetotaldensitydecreases(ρ <0).
1 1 1
hyperbolic advection equation. As a result, the nu- In addition, now (even if for small times) it is possible
merical algorithm must be designed accordingly [12]. to see the influence of the quantum corrections all over
For the equation for ψ a simple first-order upwind the sonic hole, and not just in the neighborhood of the
1
method is well suited to solve it; for the parabolic horizon. As a consequence of the shape of the quantum
equation we have implemented standard Forward-Time- sources F and G (see Fig. 3) the complex structure of
2 2
Centered-Space(FTCS)explicitmethod,aswellasstan- Fig. 4 emerges. One can see that in the region near the
dard Backward-Time-Centered-Space (BTCS) implicit horizon the fluid slows down, but there are also regions
method. Due to the short evolution time needed, the where the phonons emission induces acceleration.
limitation of the time–step required by the FTCS and
the consequent high number of iterations is not a draw-
back;inanycase,wetestedonemethodversusthe other
V. HARTLE-HAWKING STATE
andweobtainedequivalentresults. Infact,tohaveasta-
ble evolution, the time step is selected according to the
condition ∆t = α∆z2/max(ρ), since ρ is the coefficient The thermal equilibrium configuration of the Hartle-
of the ψ′′ term in the equation for ρ . In addition, for Hawking state is much simpler to treat. Since the time
1 1
the nozzle considered, we have checked that a resolution dependence drops off, the backreaction equations (23)
of∆z =2.5×10−5cm(whichcorrespondsto2000points and (24) become a simple system of algebraic equations
6
relating ρ and v : VI. CONCLUSIONS
1 1
z
A(ρ1v+ρv1)= dξF2(ξ)+const , (26) Using the continuity and Bernoulli equations, the
ZzH quantum correction (first order in ~) to a classical sta-
c2 tionary flow describing an acoustic black hole has been
A vv + ρ =G . (27)
1 ρ 1 2 evaluated in a one-dimensional approximation.
(cid:18) (cid:19)
The quantum corrections to the velocity v and to
1
The integrationconstantin Eq.(26) is chosento be zero
the density ρ profiles for the equilibrium configura-
1
inordertomakethesolutionnonsingularonthehorizon.
tion(Hartle-Hawkingstate)aredepictedinFig.5(black
The profile for v and ρ are depicted in Fig. 5 (black
1 1 lines). Thephononsbackreactioncausesthefluidtoslow
line); for the sake of comparison, we show in the same
down in the supersonic region, with the consequence of
plot the profile of v and ρ in the Unruh state for t =
1 1 a shift of the horizon to the left of the waist of the noz-
t (blue line). Inboth cases the quantumbackreaction
end zle (see Eq. (32)). In the subsonic region the velocity
correction to the velocity is positive at z = z (vertical
H increases, but the magnitude of the change is smaller
dashed line).
than the previous one. One finds a similar shape for the
In the region very close to the horizon one can make
density correction ρ , which increases in the supersonic
1
as in Ref. [6] a Taylor expansion for the background up
region and slightly decreases in the subsonic one. Fi-
to order O((z −z )5). This allows the source terms to
H nally the equilibrium temperature appears to have been
be evaluated up to linear term
lowered by the backreaction from its zero-order value
~cκ/2πκ (see Eq. (33)).
|D|A κ3 B
F2 = H −(π2+10) (28) Forthetime-dependentcase(Unruhstate)theanalysis
96π
(cid:20) hasbeenrestrictedtoveryshorttimesafterswitchingon
− π4+25π2−24κ(z−z )+O(κ2(z−z )2) , thephononsradiation. Thisbecausethequantumsource
H H
2 (F andG inEqs.(21)-(22))hasbeencomputedonlyfor
(cid:21) 2 2
A2 c2κ2 theclassicalbackground,whichjustrepresentstheinitial
G2 = H configurationoftheacousticblackhole. Amorerigorous
48π
analysis requires the time-dependence of the sources to
× (π2+6)κ(z−z )+O(κ2(z−z )2) . (29)
H H be included.
The corr(cid:2)esponding quantum corrections to th(cid:3)e velocity Within these limitations, one sees (Fig. 4) a decelera-
and to the density are tion of the fluid in the supersonic region, which causes a
drift of the horizon towards the left of the nozzle. Two
v = AHcκ2 2+π2 accelerationregionsalsoappearonbothsidesofthehori-
1 192π zon,butthe intensityofthe effect is lower. Onthe other
(cid:20)
52+35π2+π4 hand, the density correction ρ1 reflects the behaviour it
− κ(z−zH)+O(κ2(z−zH)2) , shows in the Hartle-Hawking state (on a reduced scale).
4
(cid:21)
(30)
|D|κ2
ρ = 2+π2 ACKNOWLEDGMENTS
1
192cπ
(cid:20)
+ 44−19π2−π4κ(z−z )+O(κ2(z−z )2) . ThisworkstartedwhenA.N.wasvisitingtheDepart-
H H
4 mentofAstronomyandAstrophysicsoftheUniversityof
(cid:21)
(31) Valencia. He thanks J.A. Font, L. Rezzolla and O. Zan-
ottiforsuggestionsaboutthenumericalpart;inaddition,
Setting v = v + ǫv = −c one finds the quantum
B 1 heacknowledgesthesupportofJ.Navarro–Salasandthe
corrected position of the horizon zq
H hospitality of the Department of Theoretical Physics of
theUniversityofValencia. S.F.acknowledgesthe Enrico
π2+2
zq =z − ǫA κ , (32) Fermi Center for supporting her research.
H H 192π H
which is shifted to the left of z .
H
Thequantumcorrectedequilibriumtemperaturecanalso APPENDIX A: NUMERICAL SCHEMES
be simply obtainedby evaluatingEq.(8)atz =zq with
H
v replaced by v . The result is
B In this section we report explicitly the time-evolution
~cκ ǫA algorithms. In the main text, we said that we used a
THq = 2πκ 1− 768Hπ 52+35π2+π4 κ2 , (33) standard upwind method for the equation for ψ1 and
B(cid:20) (cid:21) standardexplicitForward-Time-Centered-Space(FTCS)
(cid:0) (cid:1)
which indicates that, taking into account the backreac- or an implicit Backward-Time-Centered-Space (BTCS)
tion, the equilibrium temperature is lowered. schemes for that for ρ . In practice the upwind method
1
7
reads [13] The FTCS (explicit) and the BTCS (implicit) schemes
for the equation for ρ respectively read
1
∆t
ψn+1 =ψn −v ψn −ψn
1,i 1,i i∆z 1,i+1 1,i
(cid:0)c2 (cid:1)
+∆t − ρn +G A−1 . (A1)
ρ 1,i 2,i i
(cid:18) i (cid:19)
v ∆t ρnv′ ∆t ∆t
ρn+1 =ρn − i ρn −ρn + i i ψn −ψn −ρn ψn −2ψn +ψn
1,i 1,i 2∆z 1,i+1 1,i−1 v 2∆z 1,i+1 1,i−1 i ∆z2 1,i+1 1,i 1,i−1
i
v2v′ (cid:0) (cid:1) (cid:0) (cid:1) (cid:0) (cid:1)
+∆t − i iρn +F A−1 , (A2)
c2 1,i 2,i i
(cid:18) (cid:19)
−vi∆tρn+1 + 1+ vi2vi′∆t ρn+1+ vi∆tρn+1 =ρn +∆t F A−1+ ρnivi′ ψ1n,+i+11−ψ1n,+i−11
2∆z 1,i−1 c2 1,i 2∆z 1,i+1 1,i 2,i i v 2∆z
(cid:18) (cid:19) (cid:26) i
ψn+1 −2ψn+1+ψn+1
−ρn 1,i+1 1,i 1,i−1 . (A3)
i ∆z2
(cid:27)
where the cell index i runs from one to i . In the case selected according to the condition ∆t = α∆z2/max(ρ).
max
of the BTCS scheme ρ is obtained at every time slice For the nozzle model discussed in this paper, we have
1
(labelled by index n) as the solution of a tridiagonallin- used α = 1×10−4 to avoid stability problems and the
ear system of the form a un +b un+c un =fn that samechoicewaskeptalsofortheBTCSscheme,whichre-
i i−1 i i i i+1 i
canbeaccomplishedbyastandardlower-upper(LU)de- sultsinroughly3×104integrationsteps. Thisisnotpar-
composition of the matrix to be inverted [14]. A careful ticularlyexpensivefromthecomputationalpointofview.
treatment of the boundaries z = 0 and z = z of the Forexample,fortheresultspresentedhereweusedares-
c
numerical domain is crucial for selecting the correct so- olution of ∆z = 2.5×10−5 cm (corresponding to 2000
lution, especially when the implicit method is employed gridpoints)andthe totalevolutiontook∼11stime evo-
andsotheinversionoftheassociatedcoefficientmatrixis lution on a single-processor machine with a PentiumTM
concerned. AccordingtothephysicalmeaningoftheUn- Mprocessorat1.3GHz. Thecodewascompiledusingan
ruhstate,weimposeoutgoingconditionsatbothbound- IntelTM Fortran Compiler.
aries, i.e. un = un at i =1 and un =un at i = i We checkedconvergence of the numerical method (us-
i−1 i i+1 i max
whereun canbe either ρ orψ . Ifρ issolvedusingthe ing both BTCS or FTCS schemes) using resolutions of
1 1 1
FTCSscheme,sincethemethodisexplicitandnomatrix 500, 1000, 2000 and 4000 points, considering the case of
inversionisneeded,theproblemofsettingcorrectbound- 8000pointsasreference. Wecomputedtheerror∆f with
ary conditions is less important; in fact, it is enough to respectto the referenceresolutionas arootmeansquare
put the boundaries far enough from z to avoid any in- forf =ψ andf =ρ . Fromtherelation∆f =K∆zσ we
H 1 1
fluence on the evolution. For this kind of equations, sta- evaluatedtheconvergencerateσandweobtainedσ ≈1.3
bility has also proved to be an issue. Implementing the for both ψ and ρ . We have verified that 2000÷4000
1 1
FTCSscheme,tohaveastableevolutionthetimestepis gridpointsaresufficienttobeintheconvergenceregime.
[1] S.W. Hawking, Nature 248, 30 (1974) [7] W.G. Unruh,Phys. Rev. D14, 870 (1976)
[2] W.G. Unruh,Phys.Rev.Lett. 46, 1351 (1981) [8] J.B. Hartle and S.W. Hawking, Phys. Rev. D13, 2188
[3] R. Courant and K.O. Friedrichs Supersonic flows and (1976)
shock waves, Springer-Verlag(1948) [9] A.M. Polyakov,Phys. Lett. 103 B, 207 (1981)
[4] Artificial black holes, eds. M. Novello, M. Visser and [10] A.FabbriandJ.Navarro-Salas,Modelingblackholeevap-
G.E. Volovik, World Scientific, RiverEdge, USA(2002) oration, Imperial College Press, London (2005).
[5] C.Barcel`o, S.LiberatiandM.Visser, Analogue Gravity, [11] For a gravitational black hole the analogous expansion
gr-qc/0505065 (2005) parameterisgivenbythesquareoftheratiobetweenthe
[6] R.Balbinot,S.Fagnocchi, A.FabbriandG.P.Procopio, Plancklengthandthehorizonsize,i.e.ǫ=~/M2inunits
Phys.Rev.Lett.94,161302 (2005); R.Balbinot,S.Fag- G=c=1, see J.W. York,Phys. Rev. D31, 775 (1985)
nocchi and A.Fabbri, Phys. Rev.D71, 064019 (2005) [12] B. Gustafsson, H.O. Kreiss and J. Oliger, Time depen-
8
dent problems and difference methods, John Wiley and [14] W.H. Press, S.A. Teukolsky, W.T. Vetterling and
Sons, Inc.(1995) B.P. Flannery, Numerical Recipes, The Art of Scientific
[13] R.J. Le Veque, Numerical Methods for Conservations Computing, Cambridge UniversityPress (1992)
Laws, Birkh¨auser Verlag (1992)