Quantum simulation of Hawking radiation and curved spacetime with a superconducting on-chip black hole

Hawking radiation is one of the quantum features of a black hole that can be understood as a quantum tunneling across the event horizon of the black hole, but it is quite difficult to directly observe the Hawking radiation of an astrophysical black hole. Here, we report a fermionic lattice-model-type realization of an analogue black hole by using a chain of 10 superconducting transmon qubits with interactions mediated by 9 transmon-type tunable couplers. The quantum walks of quasi-particle in the curved spacetime reflect the gravitational effect near the black hole, resulting in the behaviour of stimulated Hawking radiation, which is verified by the state tomography measurement of all 7 qubits outside the horizon. In addition, the dynamics of entanglement in the curved spacetime is directly measured. Our results would stimulate more interests to explore the related features of black holes using the programmable superconducting processor with tunable couplers.


Introduction
In the classical picture, a particle falls into a black hole horizon and the horizon prevents the particle from turning back, then escape becomes impossible.However, taking into account quantum effects, the particle inside the black hole is doomed to gradually escape to the outside, leading to the Hawking radiation [1].The problem is that direct observation of such a quantum effect of a real black hole is difficult in astrophysics.For a black hole with solar mass, the associated Hawking temperature is only ∼ 10 −8 K and the corresponding radiation probability is astronomically small.Given by this, various analogue systems were proposed to simulate a black hole and its physical effects in laboratories [2].Over the past years, the theory of Hawking radiation has been tested in experiments based on various platforms engineered with analogue black holes, such as using shallow water waves [2][3][4][5][6][7], Bose-Einstein condensates (BEC) [8][9][10][11][12], optical metamaterials and light [13][14][15], etc.
On the other hand, the developments of superconducting processors enable us to simulate various intriguing problems of many-body systems, molecules, and to achieve quantum computational supremacy [16][17][18][19].However, constructing an analogue black hole on a superconducting chip is still a challenge, which requires wide-range tunable and site-dependent couplings between qubits to realize the curved spacetime [20].Coincidentally, a recent architectural breakthrough of tunable couplers for superconducting circuit [21], which has been exploited to implement fast and high-fidelity two-qubit gates [22][23][24][25], offers an opportunity to achieve specific coupling distribution analogous to the curved spacetime.We develop such a superconducting processor integrated with a one-dimensional (1D) array of 10 qubits with interaction couplings controlled by 9 tunable couplers, see Fig. 1, which can realize both flat and curved spacetime backgrounds.The quantum walks of quasi-particle excitations of superconducting qubits are performed to simulate the dynamics of particles in a black hole background, including dynamics of an entangled pair inside the horizon.By using multi-qubit state tomography, Hawking radiation is measured which is in agreement with theoretical prediction.This new constructed analogue black hole then facilitates further investigations of other related problems of the black hole.

Model and setup
To consider the effects of curved spacetime on quantum arXiv:2111.11092v3[quant-ph] 3 Jun 2023 matters, we consider a (1+1)-D Dirac field, of which the Dirac equation is written as (ℏ = c = 1) [26,27] iγ a e µ (a where g is the determinate of g µν , the vielbein e (a) µ satisfies the orthonormal condition e (a) µ e ν (a) = δ ν µ and the γ-matrices in the two-dimensional case are chosen to be γ = (σ z , iσ y ).In the Eddington-Finkelstein coordinates {t, x} and in the massless limit m → 0, such a Dirac field can be quantized into a discrete XY lattice model with site-dependent hopping couplings.The effective Hamiltonian reads (see Supplementary Information and ref. [20]) where σ+ j (σ − j ) is the raising (lowering) operator of the j-th qubit, µ j denotes the on-site potential, the sitedependent coupling κ j takes the form κ j ≈ f ((j − j h + 1/2)d)/4d with d being the lattice constant.Here, the function f (x) is related to spacetime metric, which is given in the Eddington-Finkelstein coordinates {t, x} as ds 2 = f (x)dt 2 − 2dtdx (see Methods and Supplementary Information).The spatial position x is discretized as x j = (j − j h )d.Since the horizon locates at f (x h ) = 0 with f ′ (x h ) > 0, the horizon in our analogues model is then defined at site j = j h where f (x h ) = 0, but the sign of κ j is different on its two sides of the horizon resulting in a black hole spacetime structure.One side of the horizon is considered as the interior of the black hole, while the opposite side represents the exterior of the black hole.
We perform the experiment to simulate the black hole using a superconducting processor with a chain of 10 qubits Q 1 -Q 10 , which represents the Hamiltonian (2), additionally with 9 tunable couplers interspersed between every two nearest-neighbour qubits, see Fig. 1.The effective hopping coupling κ j between qubits Q j and Q j+1 can be tuned arbitrarily via programming the frequency of the corresponding coupler C j , see Methods.To describe the curved spacetime experimentally, we adjust the frequencies of all the couplers, and design the effective coupling distribution as with j h = 3, ηd = 0.35 and β/(2π) ≈ 4.39 MHz.Here we choose f (x) = β tanh(ηx)/η, where η controls the scale of variation of f over each lattice site, which has the dimension of 1/d.One can verify that this function f (x) gives us a nonzero Riemannian curvature tensor and so describes a 2-dimensional curved spacetime.
As shown in Fig. 1b, the coupling κ j goes monotonically from negative to positive from Q 3 's left to right side.In this way, the information of the static curved spacetime background is encoded into the site-dependent coupling distribution.Thus, the site Q 3 where the sign of the coupling reverses can be analogous to the event horizon of the black hole, the side of negative coupling (Q 1 -Q 2 ) can be considered as the interior of the black hole, and Q 4 -Q 10 are outside the horizon.For comparison, we also realize a uniform coupling distribution with κ j /(2π) ≈ 2.94 MHz to realize a flat spacetime.In fact, from the viewpoint of the lattice qubit model, the results will be equivalent if the function κ is replaced by |κ| both in the case of curved and flat spacetime.Since we here map the coupling to the components of metric, the continuity requires κ changes the sign when passing through the analogue horizon.
In the experiment, we first prepare an initial state |ψ(0)⟩ with quasi-particle excitations, i.e., exciting qubits or creating an entangled pair.The evolution of the initial state known as quantum walk will be governed by Schrödinger equation |ψ(t)⟩ = e −i Ĥt |ψ(0)⟩ based on 1D programmable controlled Hamiltonian (2).The dynamics of the prepared states then simulate the behavior of quasi-particle in the studied (1+1)-dimensional spacetime with a designed flat or curved structure.
Quantum walks in analogue curved spacetime Figures 2a and 2b show the propagation of quasiparticles in flat and curved spacetimes, respectively.Here we initialize the system by preparing four different single-particle or two-particle states, including |ψ(0)⟩ = |1000000000⟩, |1100000000⟩, |0010000000⟩ and |0001000000⟩ with |0⟩ and |1⟩ being the eigenstates of σ+ j σ− j .Once the initial state is prepared, we apply the rectangular Z pulses on all qubits to quench them in resonance at a reference frequency of ω ref /(2π) ≈ 5.1 GHz.Meanwhile, the hopping coupling κ j between qubits is fixed as Eq.(3) (curved spacetime) or a constant (flat spacetime) by controlling couplers.After evolving for time t, all qubits are biased back to idle points for readout.The occupation of quasi-particle density distribution p j (t) := ⟨ψ(t)| σ+ j σ− j |ψ(t)⟩ is measured by averaging 5000 repeated single-shot measurements, as shown in Figs.2a and 2b.
Figure .2a shows that the propagation of quasi-particle in the flat spacetime is unimpeded, corresponding to the result of conventional quantum walk with diffusive expansion [28][29][30][31].In contrast, the particle is mainly trapped in our on-chip black hole due to the analogue gravity around the horizon Q 3 , as shown in Fig. 2b with the initial state |ψ(0)⟩ = |1000000000⟩ and |ψ(0)⟩ = |1100000000⟩.Due to the infalling Eddington-Finkelstein coordinates we took, our model only simulates the outgoing modes of the particle (see Supplementary Information).Hence, the interior and exterior of black hole are equivalent so that the same phenomenon can be observed in that case where the particle is initially prepared in the exterior of black hole (|ψ(0)⟩ = |0001000000⟩).
Here, we also present the result of the particle initialized at the horizon in Fig. 2b, i.e., |ψ(0)⟩ = |0010000000⟩.In the continuous curved spacetime, the particle initial-ized at the horizon is bound to the horizon forever due to the zero couplings on both sides of the horizon.However, in the finite-size lattice, the coupling strengths on both sides of the horizon are not strictly zero even though they are very small (≈ 0.54 MHz).The particle seems to be localized at the horizon for a very short time, but it is doomed to escape from the constraints due to the finite-size effects.
To show the accuracy of the experimental results of quantum walk in the curved spacetime, we present the fidelity F (t) = 10 j=1 p j (t)q j (t) between the measured and theoretical probability distributions p j (t) and q j (t) in Fig. 2c.The high fidelity, greater than 97% within 400 ns experiment time, implies that our experimental results are consistent with the theoretical predictions as also demonstrated by the similarity between experimental data and numerical simulations.Note that in both cases of the flat and curved spacetimes the particle will be reflected when it arrives at the boundary Q 1 or Q 10 .

Observation of analogue Hawking radiation
Black holes emit thermal radiation leading to evaporation, known as Hawking radiation.However, its observation is a challenge even for an analogue black hole due to the accuracy of the experiment.The Hawking radiation of a black hole is spontaneous in nature.The first realization of spontaneous Hawking radiation in an analogue experiment was in BEC system [8].Here we report an observation of analogue Hawking radiation on the superconducting quantum chip, which is also the first quantum realization of "lattice black hole" originally proposed by T. Jacobson more than twenty years ago [32,33].
For the initial state with a particle inside the horizon in our experiment, the evolution of the state shows the propagation of particle that results in a nonzero density of state outside the horizon is equivalent to the Hawking radiation of the black hole.Note that the Hawking radiation observed here is stimulated because we induce an excitation by flipping a qubit in |1⟩.
Defining the probability of finding a particle outside the horizon as P out = 10 j=4 p j , Fig. 2d shows a rising tendency of P out in time.This result can be considered as an important signature of Hawking radiation for the analogue black hole [3,4,14,34].
The theory of Hawking radiation points out that the probability of radiation satisfies a canonical blackbody spectrum, where E denotes the energy of particle outside the horizon, T H /(2π) = g h /(4π 2 ) is defined as the effective temperature of the Hawking radiation, and g h = 1 2 f ′ (x h ) = β/2 represents the surface gravity of the black hole [20].The derivation of Eq. ( 4) can be constructed by using the picture of quantum tunneling to obtain the tunneling rate of particle [35][36][37].We use this picture in this work to understand Hawking radiation.Such a picture is equivalent to a field theoretical viewpoint of "particle-antiparticle pairs" created around the horizon: the antiparticle (negative energy) falls into the black hole and annihilates with this particle inside the black hole, the particle outside the horizon is materialized and escapes into infinity (see Supplementary Information).Also, Eq. ( 4) can be viewed as the detailed balance relation between the creation and annihilation of particle around the horizon in a thermal environment [38,39].
The tunneling picture of Hawking radiation here is similar to the quantum fluid model of analogue horizon [40] with two differences in details.The first one is that the analogue horizon of ref. [40] is created by transonic flow but we here create analogue horizon by inhomogeneous lattice hopping.The second is that the injected beam of [40] is from the subsonic region (outside horizon) so that the reflected flow stands for the flow of Hawking radiation (classically the infalling beam should be swallowed by horizon completely and there is no reflected mode), but we here create a particle inside the horizon so the transmission flow is the Hawking radiation.
To obtain the radiation probabilities, we perform the quantum state tomography (QST) on the 7 qubits (Q 4 -Q 10 ) outside the horizon at t = 0 and t = 1000 ns, such a final time is long enough so that the particle inside the black hole has finished its tunneling to the outside but the boundary effect is negligible to the results.Here, the initial state is |ψ(0)⟩ = |1000000000⟩, i.e., a particle in the black hole has a certain position.When t = 0 ns, no radiation can be detected and all the qubits outside the horizon are almost in |0⟩, see Fig. 3a.After a long time t = 1000 ns, one may have a small chance to probe the particle outside the horizon, see Fig. 3b.The corresponding probabilities of radiation can be extracted from the measured 7-qubit density matrix.Assuming that |E n ⟩ is the n-th eigenenergy of total Hamiltonian and ρout is the density matrix outside obtained by QST, then the probability of finding a particle of energy E n outside the horizon can be obtained as P n = ⟨E n |ρ out |E n ⟩, see Methods.Although there are 2 10 = 1024 eigenstates for 10qubit Hamiltonian Eq. ( 2) and the same number of P n , the radiation states involve only 10 single-particle excited eigenstates due to the particle number conservation.As a consequence, only those P n that are corresponding to single-particle excited eigenstates have non-zero values, as shown in Fig. 3c.Therefore, we take the average of P n with the same positive energy E n as Pn to describe the average probability of finding a particle outside with E n > 0. It can be expected that the relation between Pn and E n will agree with the theoretical prediction in Eq. (4).In Fig. 3d, the simulated results show that the logarithm of the average radiation probability is approximately linear in energy with Hawking temperature 1.7 × 10 −5 K.The fitted Hawking temperature of experimental data is around ∼ 7.7 × 10 −5 K, showing validity with the same order of magnitude.The deviation between experimental data and ideal simulation data is mainly caused by the evolution of the imperfect initial state.The fidelity between the imperfect initial state in the experiment and the ideal initial state is 99.2%, which may derive from the experimental noises including XY crosstalk, thermal excitation, leakage, etc.We substitute such an experimental state for the ideal initial state in the numerical simulation of Hawking radiation, then the results of numerical simulation agree with experimental results better.
Since the analogue Hawking radiation is characterized by the temperature, we then give an estimation of how large a black hole in our real universe could reproduce the same temperature.If we consider a Schwarzschild black hole in four-dimensional spacetime with the same Hawking temperature T H , its mass can be calculated by M/M s = 6.4 × 10 −8 K/T H [1], where M s ≈ 2 × 10 30 kg is the solar mass.For the simulated black hole in our work, M/M s ∼ 10 −3 , whereas the typical value reported in BEC system for this quantity can be ∼ 10 2 [12].This significant difference in magnitude is attributed to the scales of the setup in different experimental systems.In superconducting qubits, the coupling strength is usually on the order of MHz and thus the analogue surface gravity g h is of the same magnitude, leading to T H = g h /(2π) ∼ 10 −5 K. Differently, the effective Hawking temperature of sonic black hole depends on the gradient of velocity at the analogue horizon.The BEC system and the shallow water wave system typically give us T H ∼ 10 −10 K [9-12] and 10 −12 K [4], respectively.

Dynamics of an entangled pair in the analogue black hole
Hawking predicted that the entanglement entropy increases when a black hole forms and evaporates due to the Hawking radiation.Each Hawking particle is entangled with a partner particle in the black hole.Such kind of quantum feature plays a crucial role in studying black holes and quantum information [41].
To investigate the dynamical entanglement and non-local correlation both in flat spacetime and curved spacetime, we initially prepare an entangled pair |ψ in (0)⟩ = (|00⟩ + |11⟩)/ √ 2 (Fig. 1c).The mean fidelity of prepared entangled state is up to 99.15%.The dynamics of such an initial entangled state in flat or curved spacetime is observed by time-dependent QST measurement.We obtain the two-qubit density matrix ρin (t) from the results of QST, and use it to compute the entanglement entropy and the concurrence (see Methods).In Fig. 4a, the entanglement entropy in the case of curved spacetime progressively increases due to the Hawking radiation, while in the flat spacetime it has two wavefronts resulting from the quantum interference and reflection respectively [30].On the other hand, the concurrence decreases with time in both cases, reflecting the process of entanglement being lost into the environment.However, the speed of entanglement propagation is limited by the gravitational effects near the horizon, and thus the decrease in concurrence is slowed in the curved spacetime case compared to the flat spacetime case, as shown in Fig. 4b.

Discussion
In summary, we have experimentally simulated a curved spacetime of black hole and observed an analogy of Hawking radiation in a superconducting processor with tunable couplers.An high-fidelity entangled pair is also prepared inside the horizon and the corresponding dynamics of entanglement is observed.Our results may stimulate more interests to explore the related features of black holes by using programmable superconducting processor with tunable couplers, and the techniques of calibrating and controlling coupler devices will pave the way for simulating intriguing physics with quantum many-body systems of different coupling distribution.
Our current results are a step in the direction of creating quantum systems with properties analogous to those of black holes.However, many more problems remain to be addressed in a complete simulation of quantum field theory in curved spacetime, both in theory and experiment.Theoretically, it is necessary to study different dimensional systems and investigate a comprehensive theory for mapping the various gravity fields into experimental realizable models.Experimentally, it is expected to expand the category of simulated Hamiltonians, extend the scale of qubits and enhance the control accuracy.In addition to pure analogue experiments, hybrid digital-analogue devices with substantial flexibility in near-term applications also need to be focused [42].Last but not the least, we must return to the basic problems of quantum field theory and try to translate more fundamental questions, for example, how generic is the emergence of gravity or what happens to spacetime when quantum corrections are fairly important [41].

Methods
Metric of two-dimensional spacetime.Consider a general two-dimensional spacetime background with a fixed static metric g µν , the metric in the Schwarzschild coordinates (t s , x) reads ds 2 To describe a black hole with nonzero temperature in 2-dimensional spacetime, we require that the function f has a root at x = x h with f ′ (x h ) > 0 and f (x) > 0 for x > x h standing for the exterior of the black hole, while f (x) < 0 for x < x h for the interior.The horizon of black hole then locates at x = x h .For our purpose and experimental setups, we transform above metric into "advanced Eddington-Finkelstein coordinates" by the coordinates transformation t = t s + f −1 (x)dx.In the coordinates {t, x}, the metric now becomes ds 2 = f (x)dt 2 − 2dtdx.
The differences between the "time-orthogonal coordinates" and "advanced Eddington-Finkelstein coordinates" can be found in Supplementary Note 1.
Tunable effective couplings.To construct both flat and curved spacetime background on a single superconducting quantum chip, we use tunable coupler device.
The effective coupling between nearest-neighbour qubits derives from their direct capacitive coupling and the indirect virtual exchange coupling via the coupler in between, where the former is untunable and the latter depends on the frequency of the coupler, see Supplementary Note 3. To achieve accurate control of couplings, we develop an efficient and automatic calibration for multi-qubit devices with tunable couplers, see Supplementary Note 6.In the experiments, we apply fast flux-bias Z pulses on the couplers to adjust their frequencies, contributing to the effective coupling distribution.The site-dependent coupling distribution κ j as Eq. ( 3) corresponds to the curved spacetime (β/(2π) ≈ 4.39 MHz, Fig. 1b), while a uniform coupling distribution (κ j /(2π) ≈ 2.94 MHz) is related to the flat spacetime.
Calculation of radiation probabilities.We perform the 7-qubit QST in the observation of analogue Hawking radiation and obtain the density matrix outside the horizon in the 7-qubit Hilbert space.Then we set the states of the other three qubits to |0⟩ and construct the density matrix in the 10-qubit Hilbert space ρout .The probability of finding a particle of energy E n outside the horizon P n is calculated as

Measurement of entanglement.
As shown in Fig. 1c, we prepare the initial entangled pair in the black hole by using two parallel rotations Ry π/2 and Ry −π/2 , a CZ gate ( ÛCZ = diag(1, 1, 1, −1)) and a single-qubit rotation Ry π/2 in sequence.The ideal initial state of the two qubits before the quench dynamics thus is The state of the total system (the interior of black hole and the rest) is always a pure state during the quench dynamics.Thus, the entanglement entropy of the subsystems: S(ρ in ) = S(ρ rest ), which quantifies the entanglement contained in this bipartite quantum system.In our experiment, the cost of measuring ρrest is much higher than measuring ρin due to the dimension of the Hilbert space.Therefore, we measure ρin and calculate S(ρ in ) as the entanglement measure.
To characteristic the entanglement between the two qubits in the black hole, we use the well-defined measure: concurrence [43], which can be calculated as E(ρ in ) = max{0, λ 1 − λ 2 − λ 3 − λ 4 } with λ i being the square roots of the eigenvalues of matrix ρin ρin in decreasing order, where ρin = (σ y ⊗ σy )ρ * in (σ y ⊗ σy ) is the spin-flipped state of ρin with σ y being Pauli matrix.In this paper, the spacetime geometry is present by advanced Eddington-Finkelstein coordinates (AEFC) {t, x}.Though the coordinate "t" plays the role of "time" in this system, there are a few differences compared with the usual time coordinate.In this appendix, we will give a basic introduction.
One simple way to obtain an intuition of advanced Eddington-Finkelstein coordinates {t, x} is to consider the wave propagating in flat spacetime.Let us consider a Minkowski spacetime.The usual Minkowski coordinates (MC) is {t m , x}, of which the metric reads The massless scalar field then will satisfy The solution of this equation, in general, is given by following "traveling wave solution", where ϕ 1 (t m + x) stands for the advanced solution and ϕ 2 (t m − x) stands for the outgoing solution.
To covert the advanced Eddington-Finkelstein coordinates {t, x}, we consider a coordinates transformation * These authors contributed equally to this work.† dzheng@iphy.ac.cn ‡ kaixu@iphy.ac.cn § cairg@itp.ac.cn ¶ hfan@iphy.ac.cn and so the metric then reads ds 2 = dt 2 − 2dtdx . ( It needs to note that, though dt m ̸ = dt, we still have for an arbitrary function h, i.e. the time derivatives are the same as that in the usual Minkowski coordinates and advanced Eddington-Finkelstein coordinates.Thus, the growth rate of a quantity in usual Minkowski coordinates can also be computed according to the time derivative of advanced Eddington-Finkelstein coordinates.On the contrary, if one wants to compute the spatial derivative, then the two coordinate systems will have different results in general since the derivative in left-side fixes t m but the derivative of right-side fixes t = t m + x. The propagators of wave are also very different in these two coordinate systems.From Eq. ( 3), one sees that the "traveling wave solution" in advanced Eddington-Finkelstein coordinates reads This can also be obtained from the wave equation It is a little surprising that the advanced wave ϕ 1 has no propagator!In fact the infalling mode now becomes a boundary condition rather than a propagator.If we impose a boundary condition then we have ϕ 1 (t) = 0 and there is only outgoing mode.Thus, the advanced Eddington-Finkelstein coordinates with boundary condition (9) can only represent the propagator of outgoing modes.In other words, the advanced Eddington-Finkelstein coordinates play the role of a selector to choose only outgoing modes.
Though we assume that the spacetime is flat and the matter is a scalar field in the above discussion, the basic physical picture will still be true if we consider a curved 2-dimensional spacetime and Dirac field.In Fig. 2 of our main text, one can see that our model only simulates the outgoing modes, just as we expected in the above discussion.In general, the wave in gravitational fields contains both advanced modes and outgoing modes.The Hawking radiation is an energy flux towards infinity, i.e. carried by outgoing modes.This is one reason why this paper uses advanced Eddington-Finkelstein coordinates to study Hawking radiation.

Supplementary Note 2. MORE EXPLANATION ON THE TUNNELING PICTURE OF HAWKING RADIATION
This paper uses the picture of "quantum tunneling" to understand the Hawking radiation [1][2][3].Though this is also a widespread picture to understand Hawking radiation in the community of black hole physics, it may not be familiar to other readers.Here we make a brief introduction to this picture.
At first glance, the tunneling picture is very different from the picture of "pair creation" outside the horizon.However, they are equivalent in physics.Based on the picture of "pair creation" in Hawking radiation, "particle-antiparticle pairs" can be created around the horizon.The antiparticle (negative energy) falls into the black hole and annihilates with positive energy particle inside the black hole, the particle outside the horizon is materialized and escapes into infinity.Note that the pair creation/annihilation is a virtual process, and the really materialized result is that the original particle inside the black hole disappears but an identical particle appears outside the horizon.The anti-particle of negative energy infalling the interior of the black hole can always be interpreted as a particle of positive energy outgoing from the interior, see schematic diagram Fig. S1.This leads to an equivalent picture to understand Hawking radiation via quantum tunneling: the particle inside the horizon escapes to the outside by quantum tunneling.Thus, the "tunneling picture" and "pair creation picture" are just two different pictures to understand the same physical phenomenon.Note that this "tunneling picture" does not violate causality since the spectrum is thermal and no information is carried.
Let us explain in detail how to use this tunneling picture to obtain the spectrum of radiation and corresponding temperature.For the spacetime with a black hole, the metric in the Schwarzschild coordinates {t s , x} is given by ds We consider an outgoing mode with positive energy (corresponding to the observers of infinity) of a massless scalar or Dirac field, which can be written as (ℏ = 1) By using the Eddington-Finkelstein coordinates we have Since f (x) has a root at f (x h ), we separate integration into two parts, is regular at x = x h .Here g h = f ′ (x h )/2 is the surface gravity.Thus, the outgoing positive mode then reads This outgoing mode has an infinite number of oscillations as x → x h and therefore cannot be straightforwardly extended to the inner region from the region outside the horizon.As argued in ref. [1], we can use analytic continuation to connect two branches in a complex plane: the wave function Φ describing a particle state (positive frequencies) can be analytically continued to a complex plane (see Fig. S2).Then we obtain This gives us the tunneling rate which is identical to the detailed balance relation for transition rates in a thermal environment [5].Note that the tunneling rate ( 16) stands for the rate of a single particle.It is possible that the tunneling of multiple particles happens simultaneously.For bosons we have particle number: Thus the occupation number of energy ω reads This gives us the expected distribution of bosons and the temperature reads T = g h /(2π) as predicted by Hawking.For fermions, if there is no other internal degree of freedom, the Pauli exclusion principle implies that there is at most one particle of same energy.Thus, Eq. ( 17) is replaced by particle number n : 0 1 probability: 1 P Thus the occupation number of energy ω reads This gives us the expected distribution of fermions and the temperature still reads T = g h /(2π).
In addition, the composite system of the interior of the black hole and the exterior is isolated.The interior and the exterior exchange energy and particles via the horizon.Hence, the occupations Eq. ( 18) and ( 20) can be viewed as the statistical averages of grand canonical distributions concerning bosons and fermions, respectively.can be expressed as ĤC /ℏ = where ℏ is the reduced Planck constant (for convenience ℏ will be assumed to be 1 in the following), bqj ( bcj ) and b † qj ( b † cj ) denote the annihilation and creation operators of the j-th qubit (coupler), respectively.The corresponding frequencies and anharmonicities are ω qj (ω cj ) and α qj (α cj ).Every pair of two neighbouring qubits and their middle coupler are coupled through exchangetype interactions with coupling strengths g qj ,cj , g qj+1,cj and g qj ,qj+1 .Here, the total Hamiltonian has three parts, including qubit-qubit interaction ĤQ−Q , couplercoupler interaction ĤC−C and qubit-coupler interaction ĤQ−C .The total system is equivalent to a 19-qubit Bose-Hubbard model.
In our experiment, the strong dispersive condition g qj ,cj ≪ |∆ qj ,cj | is satisfied, where ∆ qj ,cj = ω qj − ω cj is the frequency detuning.By virtue of the so-called Schrieffer-Wolff transformation one can obtain the effective qubits Hamiltonian with the corresponding dressed frequency and effective coupling strength where Λ cj qj ,qj+1 = 2/ 1/∆ qj ,cj + 1/∆ qj+1,cj is the harmonic mean of the frequencies detuning between the j-th coupler and its nearest neighbor qubits.Eq. ( 30) implies that the effective qubit-qubit coupling is derived from their direct capacitive coupling and the indirect virtual exchange coupling via the coupler in between.If the frequency of coupler is above the frequencies of qubits, Λ cj qj ,qj+1 < 0 holds so that the effective coupling g qj ,qj+1 can be tuned from positive to negative monotonically with gradually decreasing the frequency of coupler.Experimentally, we use the arbitrary waveform generator (AWG) to generate various fast-bias voltages applied to the corresponding couplers.These pulses on the Z control lines change the frequencies of couplers and then make it possible for the superconducting circuit with tunable couplers to engineer an arbitrary coupling distribution.
With g qj ,qj+1 ≪ α qj , the effective Hamiltonian Eq. ( 28) can be rewritten as a site-dependent XY model: where σ+ j (σ − j ) is the raising (lowering) operator of the j-th qubit.Here we choose For the Hamiltonian Eq. ( 31), one can map the spin variables to spinless fermion operators by introducing the Jordan-Wigner transformation [6]: ĉ † j ĉj ĉj , where the operators ĉ † and ĉj satisfy the commutation relations of fermions, i.e., {ĉ j , ĉk } = {ĉ † j , ĉ † k } = 0 and {ĉ j , ĉ † k } = δ jk .Hence, the effective Hamiltonian is mapped into a spinless fermion lattice model as By introducing a variable transformation ˆ c j (t) = (−i) j e −iµt ĉj , we obtain Here ˆ c j (t) can be viewed as a quantized operator of a discrete field φ j (t), and the spatial position can be discretized as x = x j = jd − x h , where and d denotes the lattice constant.Note the factor (−i) j is important to obtain the correct Heisenberg equation.The similar trick is widely used in artificial lattices to simulate quantum fields in flat or curved spacetimes [7][8][9].Now let us recover the continuous field φ(t, x).If we define a function f that is dependent of the spatial position x j and substitute κ j as, according to Eq. ( 35), φ(t, x j ) → ˆ c j (t)/ √ d will obey the following relation in the continuum limit, In fact, Eq. ( 37) can be considered as a special case of Dirac equation in the massless limit m → 0 if we decompose the Dirac field operator into ψ = 1 √ 2 (ξ + φ, ξ − φ) T .In the light of refs.[10,11], the Dirac equation in (1+1)dimensional curved spacetime with the metric g µν is written as where the γ-matrices in the two-dimensional case are γ 0 = σ z and γ 1 = iσ y , and the dyad is chosen as which satisfies the orthonormal condition e (a) µ e ν (a) = δ ν µ .Thus, Eq. ( 38) can be decomposed into two independent equations, In the massless limit m → 0, one can find that Eq. ( 40) is in accord with Eq. (37).Hence, what the effective Hamiltonian Eq. ( 33) describes is equivalent to a twodimensional static curved spacetime governed by the massless Dirac equation if we set κ j as There is only one single nondegenerate horizon x h so that f (x h ) = 0 and f (x) > 0 when x > x h and where g h is the surface gravity of the horizon, which gives the Hawking temperature T H = g h /(2π).In the main text, we set f (x) = β tanh ηx/η with corresponding Hawking temperature T H = β/(4π) and where η controls the scale of variation of f over each lattice site, which has the dimension of 1/d.Here, we fix j h = 3, β/(2π) = 4.39 MHz, and ηd = 0.35 in the analogue curved spacetime experiments.What we have shown above is the correspondence between XY model and the (1+1)-D Dirac field.The case of scalar field governed by Klein-Gordon equation is similar.For a complete presentation, one can refer to the earlier theoretical work [4].Here we briefly summarize the theoretical framework, as shown in Fig. S3.share one readout line equipped with a Josephson parametric amplifier (JPA) and a high-electron-mobility transistor (HEMT).Pulse on the readout transmission line is first generated as a mixture of local oscillation (LO) and the envelopes from an arbitrary waveform generator (AWG) and then demodulated by an analog digital converter (ADC).In this experiment, we replace the DC bias with a long Z square pulse generated by AWGs.Both XY and Z control signals are programmed in advance before being uploaded to AWGs.A schematic diagram of experiment setup is given in Fig. S4.
The device parameters are briefly shown in Table S1.
All the parameters are characterized by various relatively efficient and automatic methods, especially the parameters concerning couplers.Details of those methods in our experiment will be introduced in the following.

EFFICIENT AND AUTOMATIC CALIBRATION FOR MULTI-QUBIT DEVICES WITH TUNABLE COUPLERS
Before carrying out our experiment for simulating an analogue black hole, we need to calibrate all 10 qubits and find the useful parameters of 9 couplers.This is far more difficult and time-consuming than calibrating a typical 10-qubit sample without tunable couplers.In order to measure and characterize device parameters more efficiently, we adopt an automatic calibration technology based on a combination of physical models and optimization methods.

A. Spectrum of qubit and frequency calibration
First and foremost, all the qubits are individually brought up through the standard single-qubit calibration (from identifying the readout resonator frequency to calibrating π pulse).If a qubit is brought up at a certain frequency, we need to perform a two-dimensional spectroscopy measurement to extract the mapping between Z-pulse amplitude (hereinafter referred as Zpa, each unit of 200 mV) of qubit bias and its frequency (such as Fig. S5) and this will contribute to automatically calibrate all the qubits together.
According to our all-transmon sample, each transmon consists of two parallel SIS-type Josephson junctions connected by a loop which is in series with a capacitor.The critical currents of two junctions are I c1 and I c2 and E C denotes the charging energy of capacitor.By using the perturbation theory, the transition frequency can be approximately written as [12][13][14] where Φ 0 = h/(2e) is the unit flux, E JJ = I c1 I c2 Φ 0 /π denotes the geometric mean of two junctions energy in the zero field and δ = |I c1 − I c2 |/(2 I c1 I c2 ) represents the junction asymmetry.Here, the total magnetic flux Φ is in direct proportion with the strength of the magnetic field threading the loop and this weak magnetic field induced by Z pulse is approximately proportional to Zpa (Φ ∝ Zpa).Thus, the mapping between qubit Zpa and its frequency can be given by  S1.List of device parameters.Here, ω 01 q j (c j ) is |0⟩→|1⟩ transition frequency of the j-th qubit (coupler) with the corresponding readout frequency ω r q j .EC and EJJ denote the charging energy and the Josephson energy.F0,q j and F1,q j are measure fidelities of |0⟩ and |1⟩, respectively.T1,q j represents the energy relaxation time of qj at the idle point.The dephasing time T * 2,q j is characterized by the Ramsey fringe experiment, while T Echo 2,q j is measured by spin echo sequence with an inserted π pulse.The coupling strengths of exchange-type interactions between qubits and the corresponding coupler are gq j ,c j and gq j+1 ,c j , and the direct coupling of qubits is gq j ,q j+1 .FIG. S5.Experimental data of qubit automatical spectroscopy measurement.Here we take the spectrum of Q1 as an example.The black area is unscanned in order to save time.We first scan a small square area (about four columns data) around Zpa = 0 and then use polynomial curves to fit the peaks of these data.The corresponding polynomial fitting coefficients will help predict the next peak of the qubit spectrum.By constantly measuring, fitting and predicting, we obtain the experimental data of the qubit spectrum with a wide range.The mapping between qubit Zpa and its frequency can be obtained by fitting the experimental data based on Eq. (45).and Zpa(ω) ≈ arccos ± where E C can be measured by the two-photon excitation experiment (double difference between two-photon excitation frequency and qubit frequency), the remaining parameters E JJ , δ, A and ϕ will be obtained by fitting the two-dimensional spectrum of qubits.Here parameter A describes the efficiency of qubit bias which depends on the attenuation on the Z control line, while ϕ is the initial flux shift.Even if the refrigerator temperature rises and cools again, only ϕ may have some displacement.As long as the circuit wiring does not change, parameter A will keep its value.Notice, however, that Eq. ( 45) and Eq. ( 46) need to be modified by the crosstalk of Z control lines if the multi-qubit case is involved.
When we design the multi-qubit levels, Eq. ( 46) will be beneficial to obtain the corresponding Zpa according to the target frequency.A more accurate frequency calibration can be implemented by the Ramsey fringe experiment.

B. Calibration of pulse distortion and Z crosstalk
Although the Z pulse generated by AWG is designed carefully, the shape of the pulse is distorted when it interacts with the qubit.To calibrate the distortion of step response, we use several first-order infinite impulse response (IIR) filters and a finite impulse response (FIR) filter [15].Here IIR filters are designed to be an integration of several exponential functions and the FIR filter is described by a polynomial with 20 parameters.The results of pre-and post-correction are shown in Fig. S6.
For the crosstalk of Z control lines between qubits and qubits or qubits and their non-nearest neighbor couplers, a routine Z crosstalk measurement with a small scanning range is adopted, which can be used to estimate the crosstalk coefficients by measuring the frequency response to the Z control lines.However, it may be better to extend to a wider range of scanning when it comes to the Z crosstalk of couplers to qubits.If the frequency of coupler approaches the frequency of qubit, the effect of anti-crossing will be amplified due to the strong coupling between the coupler and its nearest neighbor qubit, leading to a distinctly non-linear relationship between coupler Zpa and qubit Zpa (as shown in Fig. S7b).To correct the crosstalk from classical flux crosstalk of Z control lines that basically meet the linear relationship, we first select a range of data away from the resonance points to use linear fitting, and constantly fine-tune the corresponding crosstalk coefficient until a symmetrical anti-crossing pattern is obtained.For a more accurate Z crosstalk calibration, we still take advantage of Ramsey fringe experiment, but proximity to the resonance points should be avoided.Here we emphasize that in our procedure of calibration, Z crosstalk of couplers to qubits must be corrected in order to more accurately measure the spectrum of coupler and coupling strengths, as explained in the following.

C. Spectrum of coupler and anti-crossing of energy levels
As what mentioned above, it is difficult to directly excite and measure a coupler because it has no XY control line and readout resonator.Therefore, we make use of two qubits (Q j and Q j+1 ) that are adjacent to the coupler (C j ) to perform a coupler spectroscopy measurement (after the calibration of Z crosstalk).To be specific, we apply XY excitation pulse to one of the qubits (Q j ) and vary the Zpa of coupler.If the coupler is excited to |1⟩ by the crosstalk from Q j XY line, the frequency of another qubit (Q j+1 ) will be changed due to the AC Stark effect between them.At the moment, a π pulse calibrated a, Pulse sequence for measurement of Z crosstalk of coupler to qubit.b, Before Z crosstalk is corrected, one can observe a tilted anti-crossing pattern of qubit and its nearest neighbor coupler.By constantly fine-tuning the crosstalk coefficient and applying the Zpa to compensate for crosstalk from coupler Z line, a symmetrical anti-crossing pattern will be obtained after corrected.The black area is unscanned, while the red lines are the results of the linear fitting.Here we show the experimental data of Z crosstalk calibration c, All the coefficients of Z crosstalk.Compared with the high crosstalk from couplers Z line to qubits, the absolute coefficients of Z crosstalk between qubits are all at a low level (< 2%).
before is unable to cause the perfect transition of Q j+1 due to its variation of frequency and its population of excited state will be decreased [16,17] (see pulse sequence in the inset of Fig. S8a).In this way, one can obtain the spectrum of the local Q j C j Q j+1 three-body system, which are actually the first three eigen-spectra (red lines in Fig. S8a) of the three-body Hamiltonian [16,18] with k ∈ {q j , c j , q j+1 }.In Eq. ( 47), the qubits frequencies ω qj and ω qj+1 are fixed and all anharmonicities α = −E C can be obtained by two-photon excitation measurement.Furthermore, the coupling between qubit and coupler is much stronger than g qj ,qj+1 in our device.Thus, there leaves 5 parameters to determined in this three-body Hamiltonian, namely 2 coupling strengths (i.e., g qj ,cj and g qj ,cj+1 ) and 3 parameters of ω cj (i.e., E JJ , A, ϕ in Eq. ( 45)).Note that the smallest gaps of the two anti-crossing spectral lines represent twice the coupling strengths, respectively.However, it is inaccuracy to estimate the coupling strengths only by scanning the three-body spectrum due to the broadening of spectral lines and some impure peaks [19,20].For more accurate measurement, we scan two extra anti-crossing spectrums of two-body systems Q j C j and Q j+1 C j , as shown in Fig. S8b.Truncated to two energy levels, the Hamiltonian of a qubit Q j coupled to a coupler C j can be expressed in the subspace basis {|10⟩ , |01⟩} as ĤQjCj = ω qj g qj ,cj g qj ,cj ω cj , and its eigen-energy spectra are Similarly, the eigen-energy spectra of Q j+1 C j are (50) Combining the above two equations with the diagonalization result of Eq. (47), one can finally determine the coupling strengths between coupler and qubits (i.e., g qj ,cj and g qj+1,cj ) and the mapping between coupler frequency ω cj and its Zpa.Actually, this is a multiobjective optimization problem of simultaneously fitting 3 spectroscopy results via 5 parameters.We utilize the optimization function scipy.optimize.minimize in the Python module SciPy to solve this problem.

D. Measurement of the effective coupling
To measure the effective coupling strength g qj ,qj+1 , we measure the joint probability as a function of qubitqubit swapping time t and the Zpa of coupler [19,20], as shown in Fig. S9b.Similar to Eq. ( 48), the swapping Hamiltonian of  | and coupler Zpa (or corresponding frequency) is given by each peak of normalized Fourier amplitude.The red dash line is the fitting curve of | gq j ,q j+1 | by using Eq.(54), while white dot lines denote two decoupling points ( gq j ,q j+1 = 0).As coupler frequency decreases, gq j ,q j+1 decreases from positive to zero.Once it passes the decoupling point, gq j ,q j+1 becomes negative and its absolute value will increase rapidly, especially approaching the resonance point of qubits.
which is reduced to when the two qubits are resonant, namely ω qj = ω qj+1 .Thus, the effective coupling strength can be calculated as half the Fourier frequency of probability P 01 (t).It needs to be emphasized that decoherence may cause the damping amplitude of swapping probability but does not affect the Fourier frequency.
For each Zpa of coupler (related to its frequency), we calculate g qj ,qj+1 via measuring P 01 (t) and performing Fourier transform (as shown in Fig. S9c).Subsequently, one can utilize Eq. ( 30) to draw the mapping between the effective coupling strength and coupler Zpa: where ω = ω qj = ω qj+1 is the resonant frequency of qubits, the direct coupling g qj ,qj+1 is the fitted value and the coupler frequency ω cj (Zpa) obeys Eq. (45).Hence, if the Zpa of coupler is given, the effective coupling strength can be computed via Eq.(54); or given a target coupling, one can estimate the Zpa of coupler by where E JJ and E C are the Josephson energy and the charging energy of coupler, respectively.Eq. ( 55) is a crucial foundation for engineering arbitrary coupling distribution in a superconducting circuit with tunable couplers.

Supplementary Note 7. ADDITIONAL DISCUSSION
For further discussion, we perform additional numerical simulations to compare and supplement with our results in this paper.In the following, the effects of disor-ders, different coupling distribution, finite size, and continuum limit are investigated.

A. The effects of disorders
In reality, qubits are doomed to be disturbed by various disorders, leading to the nuance between experimental conditions and theoretical assumptions.For a 1Darray of qubits, one can consider two disorders about next-nearest-neighbor (NNN) coupling g NNN and on-site potential µ with the corresponding disorder strengths W gNNN and W µ .Specifically, the Hamiltonian of disor- in the condition of strong disorder.However, we measure the NNN coupling of g j NNN ≈ 0.1 MHz and the frequencies difference of |µ j − ω ref | < 0.2 MHz with reference frequency ω ref /(2π) ≈ 5.1 GHz.According to Fig. S10, such a small degree of disorders has little impact on the results of Hawking radiation in the experiment.In fact, we measure the initial density matrix of 7 qubits outside the horizon.The fidelity between the imperfect initial state in the experiment and the ideal initial state is 99.2% (see Fig. 3a in the main text), which may be caused by the XY crosstalk, thermal excitation, leakage, etc.We substitute such an experimental state for the ideal initial state in the numerical simulation of Hawking radiation, then the results of the numerical simulation agree with the experimental results better (see Fig. 3d in the main text).

B. Different coupling distribution
Admittedly, our model does not mandate flipping the sign of coupling κ j near the horizon.In the main text, we request that the coupling goes monotonically from negative to positive (or vice versa) from the left of the horizon to its right side.This is based on the realistic consideration of the smoothness of f (x).In fact, if instead of flipping the sign of couplings they were all kept positive (or negative) inside and outside of the black hole (Fig. S11a), all of the results would be similar, as shown in Fig. S11c.For the case with non-zero coupling inside the black hole but zero coupling between all sites outside the black hole (Fig. S11b), one can find the results inside the black hole are also similar to the results in case of flipping the sign of couplings, but it is quite different for the results outside the black hole.Due to the zero coupling between all sites outside, no particles can travel in the exterior (Fig. S11d) and thus no radiation can be detected by the observer outside.

C. The finite-size effects
Here we perform the numerical simulation of a 300qubit chain to show the finite-size effects more clearly when we initialize the system by preparing a particle in the black hole.When the particle arrives at the horizon, it is going to be reflected back into the black hole in all probability but has a little chance to appear in the outside.The horizon is similar to a 'membrane' with certain transmittance, see Fig. S12.The probability of finding the particle outside P out shows a general upward trend due to the Hawking radiation.However, the particle will be reflected when it arrives at the boundary (Q 1 or Q 300 ) due to the finite-size effects.When the particle reflected by the boundary of the black hole reaches the horizon again, it has a certain probability to escape into the outside and P out thus increases again (see Fig. S12a and Fig. S12b).Conceivably, if there are no boundaries, P out will increase to a certain value and eventually the particle reaches a steady state of radiation.
In addition, the finite size can also affect the horizon.In the continuous curved spacetime, the particle initialized at the horizon is bound to the horizon forever due to the zero couplings on both sides of the horizon.However, in the finite-size lattice, the coupling strengths on both sides of the horizon are not strictly zero even though they are very small.As shown in Fig. S13a and Fig. S13b (also Fig. 2b in the main text), although the particle seems to be localized at the horizon for a very short time, it is doomed to escape from the constraints.When the particle is far from the horizon, its behavior is similar to that in flat spacetime (see Fig. S13b and Fig. S13c).

D. Continuum limit
To confirm that our model realizes an analogy of black hole, we now consider the wave prorogating in classical and continuum limit.To simulate the wave prorogating in the classical limit, we consider a state |ψ⟩ such that where ψ n = ⟨ψ|ĉ n |ψ⟩.Then Eq. ( 34) becomes We choose the function f (x) to be f (x) = tanh(αx) .
with a constant α and set µ = 0 for simplicity.This choice of metric describes an asymptotically flat spacetime since f (x) → 1 as x → ±∞.
In the asymptotically flat region x → ∞, due to the translation symmetry, we can take a planar wave ansatz  with x n = nd.Note that a coefficient 2 of n appears in Eq. ( 59), which is different from the usual mode expansion e −i(ω k t−xk+φ0) .This is because here we use Eddington-Finkelstein coordinate and need Eq. ( 59) to match the equation (11).Eq. ( 59) leads to following dispersion relationship time, we also require that the curvature radius of spacetime is much larger than the scale of wave packet so that the wave will propagate along trajectory of light.This requirement is easy found by using the Maxwell field in four-dimensional curved spacetime, Here A µ is the gauge field, □ is the d'Alembert operator in curved spacetime and R µν is the Ricci curvature tensor.In general, such an equation does not admit a plane wave solution due to the existence of the curvature term.In situations where the spacetime scale of variation of the electromagnetic field is much smaller than that of the curvature, the solutions of Maxwell's equations have a wave oscillating with nearly constant amplitude.In this case the geometric optics approximation can be used and eletromagnetic wave will travel along null geodesic approximately.
We can read the effective wavelength from Eq. ( 59 The variation of spacetime caused by the curvature in curved region of spacetime will be much lower than the variation caused by phase factors of wave.Then the geometric optics approximation can be used.In this case, the wave will travel as light rays.
In numerical simulation, we take the initial wave packet to be Gaussian with the following form Here Z is the normalized constant, ∆ describes the width of wave packet, x0 is the center of the initial wave packet and ψ n,k is defined by Eq. (59).Since the the wavepacket is modulated by Gaussian distribution, except Eq. (62) we also require ∆α ≪ 1 so that geometric optical approximation is valid.On the other hand, to have a well-defined momentum, the width of wavepacket should much be larger than the wavelength, i.e. ∆ ≫ d.
To summary, in order to check our model indeed forms an analogues of massless scalar particle in a black hole spacetime, we require the parameters to satisfy

Black Hole
Outer Space Outer Space Horizon >    >  the horizon.For finite d, the wave can only be "trapped" into the horizon for a finite time due to finite discretization.In the limit d → 0, the "light" will be trapped into the horizon forever.This is just what we expect for the black hole horizon.

E. Comparison of measured Hawking temperature with other experiments
Here we compare the measured Hawking temperature with other previous experiments, as shown in Table S2.As we mentioned in the main text, the significant difference in magnitude of the results is attributed to the scales of the setup in different experimental systems.In our superconducting qubits system, we set a stronger surface gravity, contributing to a higher Hawking temperature.

FIG. 1 .FIG. 2 .
FIG.1.On-chip analogue black hole.a, False-color image of superconducting processor and schematic analogue black hole.Ten transmon qubits, Q1 ∼ Q10, shown as crosses, are integrated along a chain with nearest-neighbour couplings.Each nearest-neighbour two qubits are coupled via a coupler, C1 ∼ C9, realized by a transmon with only a flux bias line.All the transmons are frequency-tunable, but only the qubit has the XY control line and readout resonator.The schematic image represents the background of curved spacetime simulated by this superconducting chip.The red cartoon spin located at the upper-left denotes the evolution of one quasi-particle that is initially in the black hole and the outward-going radiation.b, Schematic representation of the site-dependent effective coupling strengths κj.In the experiment, the coupling κj is designed according to Eq. (3).There is a boundary analogous to the event horizon of a black hole, where the coupling changes its sign at site Q3.Thus qubits Q1 and Q2 can be considered as the interior of the black hole, Q3 is at the horizon, and Q4 -Q10 are in the outside black hole.c, Experimental pulse sequence for observing dynamics of entanglement, which consists of three parts, i.e., (I) initialization, (II) evolution, and (III) measurement.For the initialization (I), we prepare an entangled Bell pair on Q1Q2 by combining several single-qubit pulses and a two-qubit control-phase (CZ) gate.At the left boundary of region (II), the curved (or flat) spacetime forms.Then the system will evolve according to the corresponding κj in the Hamiltonian for a time t.In region (III), we perform the state tomography measurement.

FIG. 3 .FIG. 4 .
FIG.3.Observation of analogue Hawking radiation.a, The 7-qubit density matrix at t = 0 ns.Initially, only Q1 is prepared in |1⟩ and all the qubits outside the horizon are almost in |0⟩.b, The 7-qubit density matrix at t = 1000 ns after the quench dynamics.Due to the Hawking radiation, radiation states can be detected with small probabilities.The fidelity between ideal and experimental density matrix at t = 0 and 1000 ns are 99.2% and 88.1%, respectively.c, The logarithmic probability of finding a particle outside the horizon Pn vs. its energy En.d, The logarithm of average radiation probability vs. the energy of particle when En > 0. Error bars are 1 SD calculated from the tomography data at the same energy.The slope of the red line represents the reciprocal of Hawking temperature without noise, where the Hawking temperature here is given by TH/(2π) = β/(8π 2 ) ≈ 0.35 MHz or ≈ 1.7 × 10 −5 K in Kelvin temperature.The experimental results are in agreement with the simulated data for low energy but diverge at high energy due to experiment noises.
FIG. S1.The anti-particle flow of negative energy infalling toward the interior of the black hole can always be interpreted as a particle flow of positive energy outgoing from the interior.

FIG
FIG. S2.The wavefunction are connected in the complex plane ln |x h − x| → ln |x h − x| − iπ when x runs from outside horizon into inside horizon, which yields h −x| iω/g h → |x h − x| iω/g h e πω/g h .
FIG. S4.A schematic diagram of the experimental system and partial wiring information.
FIG. S6.Experimental data of pulse shape measurement.The black area is unscanned in order to save time.Here the heatmap denotes the probability of qubit in |1⟩.The uncorrected pulse (up) is distorted, while the corrected result (down) shows a stationary step response.
FIG.S7.Experimental data of automatical Z crosstalk calibration.a, Pulse sequence for measurement of Z crosstalk of coupler to qubit.b, Before Z crosstalk is corrected, one can observe a tilted anti-crossing pattern of qubit and its nearest neighbor coupler.By constantly fine-tuning the crosstalk coefficient and applying the Zpa to compensate for crosstalk from coupler Z line, a symmetrical anti-crossing pattern will be obtained after corrected.The black area is unscanned, while the red lines are the results of the linear fitting.Here we show the experimental data of Z crosstalk calibration c, All the coefficients of Z crosstalk.Compared with the high crosstalk from couplers Z line to qubits, the absolute coefficients of Z crosstalk between qubits are all at a low level (< 2%).
FIG.S8.Experimental data of coupler automatical spectroscopy measurement.The red curves are numerical simulation results for fitting the peaks of spectroscopy data, which is based on a multi-objective optimization.a, The spectrum of the local QjCjQj+1 three-body system.The black area is unscanned, while the experimental data consists of the blue area.When the frequency of coupler is far from the qubits' frequency, we only need to scan a very narrow width like the single-qubit spectroscopy measurement.As it approaches the anti-crossing points, we increase the scan width to reduce the impact caused by predictive error, which ensures a clear three-body spectrum and saves time simultaneously.b, Experimental data of anti-crossing spectrums of QjCj (left) and Qj+1Cj (right).Here the results of C2 are taken as an example.
FIG.S9.Experimental data of the effective coupling strength measurement.a, Pulse sequence for measurement of swapping between qubits while changing the Zpa of coupler.b, Measured joint probability P01 of qubits vs Zpa of coupler (or corresponding frequency) and the swapping time.c, The Fourier transform of b, where the heatmap represents the normalized Fourier amplitude.The relation between absolute effective coupling strength | gq j ,q j+1 | and coupler Zpa (or corresponding frequency) is given by each peak of normalized Fourier amplitude.The red dash line is the fitting curve of | gq j ,q j+1 | by using Eq.(54), while white dot lines denote two decoupling points ( gq j ,q j+1 = 0).As coupler frequency decreases, gq j ,q j+1 decreases from positive to zero.Once it passes the decoupling point, gq j ,q j+1 becomes negative and its absolute value will increase rapidly, especially approaching the resonance point of qubits.
FIG. S10.The effects of two typical disorders on Hawking radiation.a, The logarithm of average radiation probability vs. the positive energy of particle with different disorder strengths of gNNN.b, The logarithm of average radiation probability vs. the positive energy of particle with different disorder strengths of µ.Here, the red solid line represents the theoretical result.
FIG.S12.Simulation of a 300-qubit chain with horizons at different locations.Here the coupling κj takes the form of Eq. (43), where d = 0.35 and β/(2π) = 4.39 MHz.From a to c, the corresponding horizons are located at Q25, Q50 and Q150, respectively.Pout is defined as the sum of probabilities of all the qubits outside the horizon.

kd ≪ 1 ,
dα ≪ 1, ∆α ≪ 1, and d/∆ ≪ 1 (64)In the Fig.S14, we set parameters ∆α = 0.2, k = 0.01, α = 0.01, d = 0.5 (a) and d = 0.05 (b) so that the geometric optical approximation is valid.From Fig.S14we see that the outgoing mode inside the black hole will be increasingly closer to the horizon but not pass through FIG. S15.Quantum walks in a 1D array of 10 superconducting qubits with black hole at the center.a, Schematic representation of the black hole at the center and the corresponding coupling distribution.b, Quantum walks in such a curved spacetime.The heatmap denotes the probabilities of excited-state for Qi in time.The horizontal axis is indexed as qubit number i, the vertical axis is time.Here we show both the numerical simulation and experiment data.c, Fidelity of the experimental data compared to ideal numerical simulation of quantum walks.

TABLE S2 .
Comparison of analogue Hawking radiation with other experiments.Considering a Schwarzschild black hole in four-dimensional spacetime with the same Hawking temperature, the mass of analogue black hole can be calculated by M/Ms = 6.4 × 10 −8 K/TH, where Ms ≈ 2 × 10 30 is the solar mass.