# 1-bit raw voltage recording system for dedicated observations of transients
at low radio frequencies

Kshitij S. Bane11affiliation: Indian Institute of Astrophysics, Koramangala
2nd Block, Bangalore 560034, Karnataka, India , Indrajit V.
Barve22affiliation: Radio Astronomy Field Station, Indian Institute of
Astrophysics, Gauribidanur 561210, Karnataka, India , G. V. S.
Gireesh22affiliation: Radio Astronomy Field Station, Indian Institute of
Astrophysics, Gauribidanur 561210, Karnataka, India , C.
Kathiravan11affiliation: Indian Institute of Astrophysics, Koramangala 2nd
Block, Bangalore 560034, Karnataka, India , and R. Ramesh11affiliation: Indian
Institute of Astrophysics, Koramangala 2nd Block, Bangalore 560034, Karnataka,
India

###### Abstract

Recently we had reported commissioning of a prototype for pulsar observations
at low radio frequencies ($<$100 MHz) using log-periodic dipole antennas
(LPDAs) in the Gauribidanur Radio Observatory
($\rm{\approx}77^{\circ}E\,14^{\circ}N$) near Bangalore in
India111https://www.iiap.res.in/?q=centers/radio. The aforementioned system
(GAuribidanur Pulsar System, GAPS) is currently being augmented to directly
digitize the radio frequency signals from the individual antennas in the
array. Our initial results using 1-bit raw voltage recording system indicates
that such a back-end receiver offers distinct advantages like, (i)
simultaneous observations of any set of desired directions in the sky with
multiple offline beams and smaller data rate/volume, (ii) archival of the
observed data with minimal resources for re-analysis in the future, either in
the same or different set of directions in the sky.

Radio observations: transients; Radio observations: pulsars; Radio
observations: instrumentation; Radio observations: low frequencies

##  1 Introduction

Latest advancements in commercial technology help to directly digitize the
radio frequency (RF) signals in radio astronomy observations and record them
as sampled voltages. Traditionally, a receiver architecture includes one or
several mixer stages performing down-conversion from RF to an intermediate or
baseband frequency, due to limited analog-to-digital (ADC) sampling rate.
However, in pace with technology development, increased ADC performance allows
for higher sampling rates. The higher the sampling rate, the higher the
frequency that can be sampled for a correct read, thereby allowing sampling at
either the observation frequency itself or a higher intermediate frequency.
Accordingly, at a sufficiently high sample rate there would be no need for
down-conversion. Eliminating the need of a mixer stage is advantageous, since
mixers in practice raise system challenges, such as introducing frequency
offsets. In addition to minimizing the number of analog components in the
receiver system and thereby the associated amplitude/phase variations, the
direct digitization offers more flexibility particularly in the efforts to
search for radio transients at low frequencies where the individual dipole
antennas have wider instantaneous sky coverage. For example, in the multi-beam
setup where beams are formed on-board digital hardware like a FPGA, the number
of beams is limited by the resources available on the FPGA. Also, once the
beams are formed in certain directions, the observations are limited to the
corresponding directions. On the contrary, if raw voltages from each antenna
are recorded, any number of beams can be formed offline to probe any desired
set of directions in the sky within the field-of-view (FoV) of the individual
antenna. Therefore, it is advantageous to have a voltage recording capability.
A high sampling rate and high resolution ADC accordingly approaches an
accurate read of an incoming signal, seemingly preferable. However, the energy
consumption at high sampling frequencies increases with the number of
quantization levels222https://github.com/bmurmann/ADC-survey. Furthermore, the
resulting data volumes and subsequent data processing/management will be
enormous in such a system. One solution to these is to record fewer number of
bits from analog-to-digital converters connected to each antenna element. Data
rates and volumes can be significantly reduced if 1-bit (sign bit) digital
receivers are employed. 1-bit digital correlators have been used in radio
astronomy (see e.g. Weinreb, 1963; Udaya Shankar & Ravi Shankar, 1990;
Nakajima et al., 1994; Ebenezer et al., 2001; Ramesh et al., 2008; Zakharenko
et al., 2016). The present work describes a 1-bit raw voltage capture system
that we have implemented in GAPS to demonstrate how such a system could be
useful for observations of pulsars and other transients at low radio
frequencies.

##  2 The antenna and analog receiver system

GAPS is a recently commissioned prototype using LPDAs for observations of the
transients in the frequency range ${\approx}$35 - 85 MHz. The antenna array
has 16 LPDAs arranged on a north-south baseline with spacing of 5 m between
the adjacent antennas (Bane et al., 2022). For the present work, we used only
8 LPDAs (see Figure 1). The half-power beam width (HPBW) of each LPDA in the
GAPS are ${\approx}80^{\circ}$ in the E-plane where the arms of the LPDA are
present, and ${\approx}110^{\circ}$ in the orthogonal H-plane. All the LPDAs
have been mounted with their H-plane in the east-west direction and E-plane in
the north-south direction. The effective collecting area and gain of each LPDA
is ${\approx}0.4{\lambda}^{2}$ (where ${\lambda}$ is the wavelength
corresponding to the observing frequency) and ${\approx}$6.5 dBi (with respect
to an isotropic radiator), respectively. The characteristic impedance of the
LPDA is ${\approx}$50 ${\Omega}$. The Voltage Standing Wave Ratio (VSWR) is
$<$ 2 over the frequency range ${\approx}$40 - 440 MHz. Since the LPDAs are
arranged on a north-south baseline, the combined response pattern of the array
(8 LPDAs) after coherent beamforming is ${\approx}110^{\circ}$ in the east-
west (hour angle/right ascenssion) direction and ${\approx}6^{\circ}$ in the
north-south (declination) direction for observations near the zenith at a
typical frequency like 60 MHz. It is a fan-beam with the above resolution in
declination. The width of the east-west response pattern is nearly independent
of frequency. Being very wide, it helps to observe a radio source continuously
for ${\approx}$7 h. This is useful since mechanical steering is difficult for
the dipole antennas used at low frequencies (see e.g. Ramesh et al., 1999b,
2012b). LMR-200 coaxial cable333https://timesmicrowave.com/wp-
content/uploads/2022/06/lmr-200-datasheet.pdf connected to the feedpoint near
the top of the LPDA transmits the RF signal incident on the LPDA to the input
of a low-pass filter with 3 dB cut-off at ${\approx}$85 MHz followed by a
high-pass filter with 3 dB cut-off at ${\approx}$35 MHz, and then a wideband
amplifier with an uniform gain ${\approx}$30 dB in the frequency range 35 - 85
MHz. The two filters and the amplifier are kept near the base of the LPDA to
minimize the length of the LMR-200 cable and hence the transmission loss. The
high- and low-pass filters help to attenuate the unwanted signal at
frequencies ${\lesssim}$30 MHz and ${\gtrsim}$85 MHz. Both the filters have
steep reduction in gain beyond their cut-off frequencies. For e.g. the
amplitude of the input signal is reduced by ${\approx}$30 dB at 25 MHz and 95
MHz in the aforementioned high- and low-pass filter, respectively. Due to
this, amplification of any radio frequency interference (RFI) present outside
the frequency range of our observation (i.e. 35 - 85 MHz) is very minimal. The
passband is relatively free of RFI (see e.g. Bane et al., 2022). The output
corresponding to each LPDA are then transmitted to a central cabin near the
center of the array where they are again high-pass filtered and amplified.
Note that the 2nd high-pass filter in the signal chain further suppresses the
RFI at frequencies ${\lesssim}$30 MHz. The RF signals are then transferred to
the receiver building (located ${\approx}$300 m away) independently via eight
separate, low-loss LMR-400 coaxial cables444https://timesmicrowave.com/wp-
content/uploads/2022/06/lmr-400-datasheet.pdf of identical RF characteristics.
The lengths of all the cables are the same and they were also phase equalized.
The cables are buried ${\approx}$1 m below the ground level to minimize
possible diurnal variations in their characteristics.

Figure 1: GAPS array configuration used in the present work.

The array was initially equipped with an analog beamformer network (Bane et
al., 2022) where phase/delay cables controlled by diode switches were used to
observe a particular direction of the sky at any given epoch (see e.g.
Landecker, 1984; Ramesh et al., 1998). Different lengths of phase/delay cables
must be used in the latter to observe different regions of the sky. This
introduces systematic gain variations. Furthermore, simultaneous observations
are practically difficult. Once the beam is formed in a particular direction,
observations will be restricted to that direction. These are common to any
other similar analog hardware used elsewhere. Note that we used only 8 LPDAs
in the GAPS for the present work since our aim is to understand the digital
beamformer described in the present work for future use with the Gauribidanur
RAdioheliograPH (GRAPH) which has 8 LPDAs per group (Ramesh, 2011; Ramesh et
al., 2014), and use GRAPH for dedicated observations of non-solar transients
particularly during the local night time. Further, the analog-to-digital
convertors available with us can also handle 8 inputs only. Note that the
GRAPH is a T-shaped radio interferometer array consisting of 512 LPDAs. The
individual arms of the array are oriented along the east-west and south
directions. The effective collecting area ($\rm A_{e}$) of the array is
$\rm{\approx}6400m^{2}$ at a typical frequency like 60 MHz. The sky coverage
is ${\approx}$8800 square degrees, independent of frequency. Both these
numbers are reasonably larger, an important requisite for observations of non-
solar transients. For comparison, the fully completed NenuFAR in France is
expected to have $\rm A_{e}{\approx}20000m^{2}$ and sky coverage
${\approx}100$ square degrees, at 60 MHz555https://nenufar.obs-
nancay.fr/en/astronomer/. The corresponding numbers for LOFAR which is
currently operational are $\rm A_{e}{\approx}9600m^{2}$ and sky coverage
${\approx}19600$ square degrees, at 60 MHz (van Haarlem et al., 2013). The
aforementioned characteristics of GRAPH, along with the prospects of dedicated
and continuous observing period of ${\approx}$14 h every day exclusively for
non-solar transients suggests that GRAPH could be an useful instrument.

##  3 The digital receiver system

Figure 2 shows the schematic of the 1-bit raw voltage recording receiver. The
digital receiver is implemented on Reconfigurable Open Architecture Computing
Hardware (ROACH) from the Collaboration for Astronomy Signal Processing and
Electronics Research (CASPER). The ROACH board hardware has the Xilinx
Virtex-5 FPGA666https://casper.astro.berkeley.edu/wiki/ROACH. The RF signals
from the each of the eight antennas pass through independent amplifier and
bandpass filter of 50 - 70 MHz before being fed to ADCs. We have used two
Quad-ADC cards, each containing four AD9480
ICs777https://www.analog.com/en/products/AD9480.html. The signals are sampled
at 90 MHz. A single clock is used to sample all the antennas. The sampled band
lies in the second Nyquist zone of 90 MHz. The ADC converts the input voltage
to a 8-bit Fixed point number (Fix 8.7) between -1 and +1. The sign bit from
each of these numbers is separated. So we have 1-bit from each antenna
channel. Therefore, the ADC output is an unsigned 8-bit number (UFix8.0)
representing the input voltage from each of the 8 channels. Such 1024x8 8-bit
numbers are packetized together and sent to the recording computer over 10
Gigabit Ethernet (10GbE). Note that 10GbE allows transmission of 64-bit
numbers888https://casper-
toolflow.readthedocs.io/projects/tutorials/en/latest/tutorials/roach/tut{_}ten{_}gbe.html.
So, one packet has 1024 64-bit numbers. This results in a data rate of
${\approx}$0.72 Gbits/s. This is much smaller compared to the the data rate of
${\approx}$5.6 Gbits/s which would have resulted if 8-bits were recorded from
each antenna channel. To avoid packet losses while recording, the data is
captured using n2disk program999https://www.ntop.org/products/traffic-
recording-replay/n2disk/ on a 32 GB RAM computer. Each packet also contains a
user datagram protocol (UDP) header(58 bytes) and a custom header (80 bytes)
that contain the observation details, 1 pulse per second (PPS) count, and a
packet count. The 1 PPS signal is derived from a GPS
clock101010https://timing.trimble.com/wp-content/uploads/thunderbolt-e-gps-
disciplined-clock-datasheet.pdf and given to the FPGA via the ADC card to
generate the 1 PPS count. The packet-count is a unique packet number assigned
to each packet. It is used to check for any packet losses. The recorded file
contains a time series of 1-bit voltages from each antenna channel. Some of
the characteristics are listed in Table 1.

Figure 2: 1-bit raw voltage recording system used in the GAPS. Table 1: Characteristics of GAPS Number of channels | 8
---|---
Sampling rate | 90 MSPS
Sampled bits | 1
Observation band | 45 - 90 MHz (50 - 70 MHz for the present work)
Beamforming | offline
FFT size | Variable (8192 for the present work)
Frequency resolution | 10.98 kHz
Post-detection time resolution | Variable (min 91.02 ${\mu}$s for one spectrum)
Declination coverage |  $-26.4^{\circ}$ S to $+53.6^{\circ}$ N
1 | 2

##  4 Offline beam formation

Figure 3 shows the offline beamforming pipeline. First, the packets are
unscrambled to get voltage time series of the individual antenna channels.
Delays are applied to these time series according to the desired declination
where the beam needs to be formed. Note that the array is oriented in the
north-south direction, so the beam forming is performed along the declination.
Based on the declination, the geometric delays are calculated using the
formula $t_{j}$=$d_{j}sin{\theta}/c$, where $t_{j}$ is the delay required for
$j^{th}$ antenna at a distance of $d_{j}$ from the reference antenna for the
source at angle ${\theta}$, and $c$ is the velocity of light. The sampling
interval ($t_{s}$) is 1/90 MHz = 11.11 ns. So, the time series are delayed by
the closest integer multiple of 11.11 ns (i.e. $mt_{s}$, where $m$ is an
integer). The uncompensated fractional delays (${\Delta}t_{j}$=$t_{j}$ \-
$mt_{s}$) can be compensated by applying a phase gradient in the spectral
domain, i.e. ${\phi}_{j,k}$=$2{\pi}F_{k}{\Delta}t_{j}$, where $F_{k}$ is
$k^{th}$ frequency bin. The delayed voltage streams are Fourier transformed
using FFT. Subsequenty, the phase compensation is performed by multiplying the
complex Fourier output by $e^{-i{\phi_{j}}}$. An 8192-point FFT is performed,
resulting in 4096 positive frequency bins with a resolution of 10.986 kHz. The
number of FFT points and hence the resolution can be changed depending on the
requirement since the beamforming is performed offline. After the phase
compensation, the individual spectra from each antenna are added together,
modulo-squared, and integrated according to the desired temporal resolution.
Multiple beams pointing at different declinations can be formed by changing
the delay and phase values. A ‘.SPEC’ file is stored for each beam containing
the corresponding power spectrum series. These are further processed using the
standard PulsaR Exploration and Search TOolkit (PRESTO, Ransom, 2011).

Figure 3: Offline beamforming pipeline used in the GAPS.

##  5 Observations

Figure 4: Observations of the pulsar B0834+06 (J0837+0610) on 2023 October 21
with GAPS in the frequency range of 50 to 70 MHz for a period of ${\approx}$4
h. The data were analyzed using PRESTO.

Figure 5: Average profiles of the pulsars observed with GAPS using 1-bit raw
voltage recording system and offline beamformation.

The initial run of the new system was carried out on 2023 August 10 for
${\approx}$5 h, starting from ${\approx}$08:49 LST. Three beams pointing
towards declinations $6.17^{\circ}$ N, $7.92^{\circ}$ N, and $15.85^{\circ}$ N
were formed offline from the raw voltage data recorded. The power spectrum
series from these three beams were stored in three separate .SPEC files with
time resolutions of 4 ms, 2 ms and 1 ms, respectively. The files were
converted to SIGPROC ‘filterbank’ (.fil)
format111111http://sigproc.sourceforge.net/sigproc.pdf. As first step, RFI
excision was performed using rficlean (Maan et al., 2021) to mitigate strong
periodic RFI. Subsequently we used rfifind from PRESTO, which searches for
time and frequency domain RFI and creates mask files. This is followed by
prepfold which is used to fold and dedisperse the data and detect pulsars.
Three pulsars B0834+06 (J0837+0610), B0950+08 (J0953+0755) and B1133+16
(J1136+1551) were detected in the above three beams, respectively (see Figures
4 & 5). A solar event (type III burst, McLean & Labrum, 1985; Ramesh &
Ebenezer, 2001) was also recorded during the observation at ${\approx}$11:03
LST. The declination of Sun was ${\approx}15.7^{\circ}$ N. A separate beam was
formed pointing towards the above declination with time resolution of
${\approx}$10 ms. We carried out analysis with single pulse search pipeline
developed using PRESTO routines. The maximum SNR was obtained for dispersion
measured (DM) ${\approx}$ 13.2 $\rm pc\,cm^{-3}$. The de-dispersed dynamic
spectra of the solar burst is shown in Figure 6. The burst was noticed in the
observations carried out with the solar radio spectrographs in the
Gauribidanur observatory also (Benz et al., 2009; Kishore et al., 2014).

To verify the performance of GAPS, we formed 41 beams covering the declination
range from $-26.3^{\circ}$ S to $+53.7^{\circ}$ N in steps of $2^{\circ}$ for
the duration of the solar burst in Figure 6. Mean count during the solar burst
at 60 MHz in each of the 41 beams were obtained. The count is maximum for the
beam closest to the Sun’s declination. A Gaussian fit to the data points
indicates HPBW ${\approx}6^{\circ}$ (see Figure 7), which matches with the
predicted HPBW of the GAPS antenna array in the present case (see Section 2).
On 2023 October 18, a 4 h observation run starting at ${\approx}$17:27 LST was
carried out. A beam was formed pointing towards declination
${\approx}21.88^{\circ}$ N with 4 ms time resolution. Pulsar B1919+21
(J1921+2153) was detected in this data (see Figure 5). A subsequent
observation for a period of ${\approx}$7.2 h was performed on 2023 October 21
starting at ${\approx}$06:36 LST. Pulsars B0834+06 (J0837+0610), B0919+06
(J0922+0638), B0950+08 (J0953+0755), B0943+10 (J0946+0951) and B1133+16
(J1136+1551) were detected using ${\approx}$6 h data around the transit of
each pulsar (see Figure 5). Four beams were formed pointing towards
declinations $6.17^{\circ}$ N, $7.92^{\circ}$ N, $9.86^{\circ}$ N and
$15.85^{\circ}$ N with 1 ms resolution. Pulsars B0834+06 and B0919+06 were
detected in the same beam. The estimated parameters of the different pulsars
mentioned above are listed in Table 2. The corresponding numbers reported
previously by Bondonneau et al. (2020) from online, multi-bit digitization
observations elsewhere in nearly the same frequency range are mentioned within
the brackets. The period & DM for each pulsar in the two cases agree well. The
pulsar average profiles in Figure 5 also closely match with the average
profiles at low frequencies reported in the literature (Malov & Malofeev,
2010; Hassall et al., 2012; Zakharenko et al., 2013; Stovall et al., 2015;
Pilia et al., 2016; Bilous et al., 2020; Bondonneau et al., 2020). Any minor
differences are likely due to noisy detections in the present case. For
example in the case of Bondonneau et al. (2020), the bandwidth and $\rm A_{e}$
used were higher than the present observations by factors of ${\approx}$2.75 &
${\approx}$10, respectively. Furthermore, 12-bit ADC was used in their
observations. Its quantization efficiency is ${\approx}$(1/0.64) times higher
than that of the 1-bit ADC used in the present observations (Thompson et al.,
2001). The duration of observation for the respective pulsars in the two cases
are also different. If we take into consideration all these factors, then the
expected SNR for the pulsar B0834+06 in GAPS observations as compared to
Bondonneau et al. (2020) is,
${\frac{309*0.64}{10{\sqrt{2.75/6}}}}{\approx}30$. Compared to this, the
observed SNR with GAPS is ${\approx}$36\. Note that the duty cycle mentioned
in Table 1 corresponds to the effective pulse width in pulse profiles at 50%
of the highest peak intensity (i.e. w50).

Figure 6: De-dispersed dynamic spectrum (50 to 70 MHz) corresponding to the
type III solar radio burst observed with the GAPS on 2023 August 10. The DM
used was $\rm{\approx}13.2\,pc\,cm^{-3}$. The time in the x-axis starts from
${\approx}$11:02:30 LST, and the interval is 10 s. The upper panel show the
spectrally averaged time profile of the burst. The right panel shows the
temporally averaged spectral profile of the burst. The dotted lines in the
upper panel indicate the main phase of the burst.

Figure 7: The declination beam of GAPS at 60 MHz obtained by using multiple
beams formed during the solar observations in Figure 6.

##  6 Conclusions

We have presented a 1-bit raw voltage recording system for observations of
pulsars at low radio frequencies. Other than reduction in the senstivity which
could be partly compensated by the use of higher sampling rate (Van Vleck &
Middleton, 1966; Burns & Yao, 1969; Thompson et al., 2001), the 1-bit
quantization does not have an effect on the parameters of a pulsar like
period, DM, etc. (see Table 2). Most of the modern low-frequency radio
telescopes have multi-bit digital receivers which provide higher dynamic range
(Taylor et al., 2012; van Haarlem et al., 2013; Prabu et al., 2015; Reddy et
al., 2017; Zarka et al., 2020; Girish et al., 2023). However, in a relatively
minimal RFI location like Gauribidanur (Monstein et al., 2007; Kishore et al.,
2015; Hariharan et al., 2016b; Bane et al., 2024), a lower dynamic range
digital receiver can be used (see e.g. Zakharenko et al., 2016). Our results
indicate that in such circumstances the sensitivity to observe faint events
within a short time interval depends largely on the effective collecting area
and the observing bandwidth. It is well known that absolute flux measurments
can be difficult in 1-bit receivers (Van Vleck & Middleton, 1966; Udaya
Shankar & Ravi Shankar, 1990; Ramkumar et al., 1994; Ramesh et al., 2006b;
Stein, 2019). But in a wide-field multi-beam system where the primary
requirement is to cover maximum area of the sky with high temporal and
spectral resolution, such binary receiver systems can be beneficial. Overall,
the 1-bit receiver offers flexibility, allowing multiple simultaneous beams in
different directions, different integration times and frequency channel
widths, and different filtering algorithms to be used. That wouldn’t make your
purposes less compelling but it would certainly advocate for this kind of
recording mostly being used for bright events. Further, the archival data can
be stored in raw form with fewer resources and reprocessed in the future as
per requirements. Reports indicate that 1-bit raw voltage recording had been
used at higher frequencies also (Teng et al., 2015). With the currently
available technology, sampling rates upto ${\approx}$40 GHz can be achieved.

Table 2: Characteristics of the pulsars observed with GAPS using 1-bit raw voltage recording system (Numbers within the brackets are from Bondonneau et al. (2020)) Pulsar | Period | DM | Duty cycle | SNR | Observing
---|---|---|---|---|---
| (s) | ($\rm pc/cm^{3}$) | (%) | (folded) | duration
|  |  |  |  | (h)
B0834+06 | 1.274${\pm}0.000043$ | 12.873${\pm}0.05$ | 6.4${\pm}0.8$ | 36.5 | 6
| (1.274) | (12.864) | (6.9) | (309) | (1)
B0919+06 | 0.431${\pm}0.000016$ | 27.330${\pm}0.02$ | 14.4${\pm}0.8$ | 6.0 | 6
| (0.431) | (27.296) | (15.4) | (144) | (3)
B0943+10 | 1.098${\pm}0.000069$ | 15.348${\pm}0.02$ | 12.9${\pm}0.8$ | 6.3 | 6
| (1.098) | (15.329) | (15.2) | (148) | (2.5)
B0950+08 | 0.253${\pm}0.000050$ | 2.965${\pm}0.05$ | 15.3${\pm}0.8$ | 24.6 | 6
| (0.253) | (2.971) | (14.6) | (140) | (1)
B1133+16 | 1.188${\pm}0.000023$ | 4.848${\pm}0.03$ | 12.7${\pm}0.8$ | 8.4 | 6
| (1.188) | (4.846) | (18.7) | (261) | (2)
B1919+21 | 1.337${\pm}0.000017$ | 12.420${\pm}0.02$ | 7.8${\pm}0.8$ | 12.5 | 4
| (1.337) | (12.437) | (8.4) | (180) | (1)

We are grateful to Gauribidanur Observatory team for their help in the
observations and upkeep of the facilities. A. A. Deshpande is thanked for his
encouragement and suggestions. We acknowledge the referee for his/her kind
comments, which helped to present the results more clearly.

## References

  * Bane et al. (2022) Bane, K. S., Barve, I. V., Gireesh, G. V. S., Kathiravan, C., & Ramesh, R. 2022, J. Astron. Teles. Instr. Sys., 8, 017001
  * Bane et al. (2024) —. 2024, J. Astron. Teles. Instr. Sys., 10, 014001
  * Benz et al. (2009) Benz, A. O., Monstein, C., Meyer, H., et al. 2009, Earth, Moon and Planets, 104, 277
  * Bilous et al. (2020) Bilous, A. V., Bondonneau, L., Kondratiev, V. I., et al. 2020, Astron. Astrophys., 635, A75
  * Bondonneau et al. (2020) Bondonneau, L., Grießmeier, J.-M., Theureau, G., et al. 2020, Astron. Astrophys., 635, A76
  * Burns & Yao (1969) Burns, W. R., & Yao, S. S. 1969, Radio Sci., 4, 431
  * Ebenezer et al. (2001) Ebenezer, E., Ramesh, R., Subramanian, K. R., Sundara Rajan, M. S., & Sastry, C. V. 2001, A&A, 367, 1112
  * Girish et al. (2023) Girish, B. S., Reddy, S. H., Sethi, S., et al. 2023, J. Astrophys. Astron., 44, 1641011
  * Hariharan et al. (2016b) Hariharan, K., Ramesh, R., Kathiravan, C., Abhilash, H. N., & Rajalingam, M. 2016b, ApJS, 222, 21
  * Hassall et al. (2012) Hassall, T. E., Stappers, B. W., Hessels, J. W. T., et al. 2012, Astron. Astrophys., 543, A66
  * Kishore et al. (2014) Kishore, P., Kathiravan, C., Ramesh, R., Rajalingam, M., & Indrajit, V. B. 2014, Sol. Phys., 289, 3995
  * Kishore et al. (2015) Kishore, P., Ramesh, R., Kathiravan, C., & Rajalingam, M. 2015, Sol. Phys., 290, 2409
  * Landecker (1984) Landecker, T. L. 1984, IEEE Trans. Instr. Meas., 78
  * Maan et al. (2021) Maan, Y., van Leeuwen, J., & Vohl, D. 2021, Astron. Astrophys., 650, A80
  * Malov & Malofeev (2010) Malov, O. I., & Malofeev, V. M. 2010, Astron. Rep., 54, 210
  * McLean & Labrum (1985) McLean, D. J., & Labrum, N. R. 1985, in Solar Radio Physics: Studies of emission from the Sun at metre wavelengths, Cambridge University Press
  * Monstein et al. (2007) Monstein, C., Ramesh, R., & Kathiravan, C. 2007, Bull. Astron. Soc. India, 35, 473
  * Nakajima et al. (1994) Nakajima, H., Nishio, M., Enome, S., et al. 1994, Proc. IEEE., 82, 705
  * Pilia et al. (2016) Pilia, M., Hessels, J. W. T., Stappers, B. W., et al. 2016, Astron. Astrophys., 586, A92
  * Prabu et al. (2015) Prabu, T., Srivani, K. S., Roshi, A. D., et al. 2015, Exp. Astron., 39, 73
  * Ramesh (2011) Ramesh, R. 2011, in 1st Asia-Pacific Sol. Phys. Meeting, ed. A. R. Choudhuri & D. Banerjee, Astron. Soc. India Conf. Ser. 2, 55
  * Ramesh & Ebenezer (2001) Ramesh, R., & Ebenezer, E. 2001, ApJ, 558, L141
  * Ramesh et al. (2012b) Ramesh, R., Kathiravan, C., Indrajit, V. B., & Rajalingam, M. 2012b, ApJ, 744, 165
  * Ramesh et al. (2014) Ramesh, R., Kathiravan, C., Sundara Rajan, M. S., Indrajit, V. B., & Rajalingam, M. 2014, in The Metrewavelength Sky, ed. J. N. Chengalur & Y. Gupta, Astron. Soc. India Conf. Ser. 13, 19
  * Ramesh et al. (2008) Ramesh, R., Kathiravan, C., Sundara Rajan, M. S., Indrajit, V. B., & Sastry, C. V. 2008, Sol. Phys., 253, 319
  * Ramesh et al. (1999b) Ramesh, R., Subramanian, K. R., & Sastry, C. V. 1999b, Sol. Phys., 185, 77
  * Ramesh et al. (1998) Ramesh, R., Subramanian, K. R., Sundara Rajan, M. S., & Sastry, C. V. 1998, Sol. Phys., 181, 439
  * Ramesh et al. (2006b) Ramesh, R., Sundara Rajan, M. S., & Sastry, C. V. 2006b, Exp. Astron., 21, 31
  * Ramkumar et al. (1994) Ramkumar, P. S., Prabu, T.and Girimaji, M., & Markendeyalu, G. 1994, J. Astrophys. Astron., 15, 343
  * Ransom (2011) Ransom, S. 2011, Astrophysics Source Code Library, record ascl:1107.017
  * Reddy et al. (2017) Reddy, S. H., Kudale, S., Gokhale, U., et al. 2017, J. Astron. Instr., 6, 1641011
  * Stein (2019) Stein, M. S. 2019, arXiv:1905.12528
  * Stovall et al. (2015) Stovall, K., Ray, P. S., Blythe, J., et al. 2015, Astrophys. J., 808, 156
  * Taylor et al. (2012) Taylor, G. B., Ellingson, S. W., Kassim, N. E., et al. 2012, J. Astron. Instr., 1, 1250004
  * Teng et al. (2015) Teng, H. F., Zhang, U. H., Chiueh, T. H., et al. 2015, IEEE Trans. Instr. Meas., 64, 299
  * Thompson et al. (2001) Thompson, A. R., Moran, J. M., & Swenson Jr., G. W. 2001, Interferometry and Synthesis in Radio Astronomy 2nd Edition (NewYork: Wiley-Interscience)
  * Udaya Shankar & Ravi Shankar (1990) Udaya Shankar, N., & Ravi Shankar, T. S. 1990, J. Astrophys. Astron., 11, 297
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, Astron. Astrophys., 556, A2
  * Van Vleck & Middleton (1966) Van Vleck, J. H., & Middleton, D. 1966, Proc. IEEE, 54, 2
  * Weinreb (1963) Weinreb, S. 1963, Ph.D Thesis, MIT Tech. Rep. 412
  * Zakharenko et al. (2016) Zakharenko, V., Konavalenko, A. A., Zarka, P., et al. 2016, J. Astron. Instr., 5, 1641010
  * Zakharenko et al. (2013) Zakharenko, V. V., Vasylieva, I. Y., Konovalenko, A. A., et al. 2013, Mon. Not. Roy. Astron. Soc., 431, 3624
  * Zarka et al. (2020) Zarka, P., Denis, L., Tagger, M., et al. 2020, in URSI GASS 2020, Session J01 New Telescopes on the Frontier, https://tinyurl.com/ycocd5ly


# A quantum walk simulation of extra dimensions with warped geometry

Andreu Anglés-Castillo  Universitat de València-CSIC, Departament de Física
Teòrica & IFIC, Burjassot (València), 46100, Spain<EMAIL_ADDRESS>Armando Pérez  Universitat de València-CSIC, Departament de Física Teòrica &
IFIC, Burjassot (València), 46100, Spain

###### Abstract

We investigate the properties of a quantum walk which can simulate the
behavior of a spin $1/2$ particle in a model with an ordinary spatial
dimension, and one extra dimension with warped geometry between two branes.
Such a setup constitutes a $1+1$ dimensional version of the Randall-Sundrum
model, which plays an important role in high energy physics. In the continuum
spacetime limit, the quantum walk reproduces the Dirac equation corresponding
to the model, which allows to anticipate some of the properties that can be
reproduced by the quantum walk. In particular, we observe that the probability
distribution becomes, at large time steps, concentrated near the "low energy"
brane, and can be approximated as the lowest eigenstate of the continuum
Hamiltonian that is compatible with the symmetries of the model. In this way,
we obtain a localization effect whose strength is controlled by a warp
coefficient. In other words, here localization arises from the geometry of the
model, at variance with the usual effect that is originated from random
irregularities, as in Anderson localization. In summary, we establish an
interesting correspondence between a high energy physics model and
localization in quantum walks.

## Introduction

Quantum walks (QWs) constitute an interesting possibility for simulating
physical phenomena from many fields. The discrete time version describes the
motion of a spin $1/2$ particle on a lattice. For instance, by simply
incorporating suitable position-dependent phases on the unitary operator that
implements the time evolution, one can mimic the effects of an external
electromagnetic field [1, 2, 3, 4, 5, 6, 7, 8]. In the continuum limit (when
both the time step and the lattice spacing tend to zero), the Dirac equation
in presence of such fields is recovered. In an analogous way, the motion of a
Dirac particle in presence of a gravitational field can be simulated by an
appropriate choice of the operator that drives the evolution, either on a
rectangular or other types of lattices [3, 9, 10]. Other scenarios include
vacuum or matter neutrino oscillations [11, 12, 13], and one can even
establish some connections to lattice field theories [14].

There is also a different connection of QWs with quantum field theories,
namely the possibility to explore some models which include extra dimensions,
which are only manifested at very high energies. The possibility of extra
dimensions of space was first suggested by Theodor Kaluza and Oscar Klein [15,
16] seeking an unified theory of electromagnetic and gravitational fields into
a higher dimensional field, with one of the dimensions compactified.
Experimental data from particle colliders restrict the compactification radius
to such small scales that it becomes virtually impossible to explore these
extra dimensions. Different ideas have been proposed to overcome this
difficulty, for example the domain wall model introduced by Rubakov and
Shaposhnikov [17], in which the particle couples to an external scalar field.
The motion of a spin $1/2$ particle moving inside such a geometry was analyzed
in [18]. In addition to recovering the corresponding Dirac equation in the
continuum limit, the QW shows, at finite spacetime spacing, localization of
the particle within the brane due to the coupling to the field.

Spatial localization is an important phenomenon in physics, which appears
within the context of diffusion processes in lattices. It can arise from
random noise on the lattice sites, giving rise to Anderson localization [19]
and causing a metal-insulator transition, but it can also be the consequence
of the action of an external periodic potential (see e.g. [20, 21, 22]).
Similarly, one obtains localization for the 1-dimensional QW when spatial
disorder is included [23, 24, 25], non-linear effects [26], or by the use of a
spatially periodic coin [27]. The results in [18] show, however, that
localization can also appear as a consequence of the interaction with a smooth
external potential, instead of a random, or even periodic, perturbation.

In this paper, we investigate localization effects that arise within a
different context, which is also inspired on high energy physics, and was
originally proposed to address the hierarchy problem (the observed difference
between the Higgs mass, and the Planck scale, in many orders of magnitude),
and is commonly referred to as the Randall-Sundrum model [28]. This model
assumes an extra dimension which extends between two branes (with a topology
that will be discussed later). Here we consider a simplified version with one
ordinary spatial dimension and one extra dimension, and define a QW that
reproduces the dynamics of a spin $1/2$ particle in the continuum spacetime
limit.

Unlike the Rubakov and Shaposhnikov model, there is no coupling to an external
scalar field. Instead, this model presents a warped geometry along the extra
dimension. As we will show, this curvature is at the root of a localization
effect of the QW towards the second (low energy) brane. The stationary states
of the model in the continuum limit become concentrated close to the low
energy brane for high values of the warp coefficient, which quantifies the
strength of the localization. The localization of the QW can be analyzed by
quantifying its overlap with these stationary states. This allows us to tailor
the dynamics of the QW, showing a different behavior as the value of the warp
coefficient is changed. In this way, we arrive at a QW model with a rich
phenomenology, where some properties are inherited from the continuum field
theoretic model. There is, in this sense, a mutual multidisciplinary benefit:
one can design a QW which simulates an important high energy physics model. In
exchange, the knowledge of the continuum properties is useful to understand,
and to control, the dynamics of the QW in different regimes.

This paper is organized as follows. We first define the Randall-Sundrum model
in $1+1$ spatial dimensions, along with its main properties. We pay special
attention to the stationary states of the Hamiltonian, which play a crucial
role in understanding the dynamics of the proposed QW. Next, we define a QW
which allows to recover the dynamics of the Randall-Sundrum model for a spin
$1/2$ particle, and we study its phenomenology. Namely, we show that the
distribution probability, as well as the expected value of the position along
the extra dimension, approaches the lower brane at large time, and that this
approaching proceeds more slowly for larger values of the warp coefficient,
which turns out to be the main parameter in controlling the dynamics. We also
analyze the entanglement entropy between spatial and internal degrees of
freedom, exhibiting a complex behavior as a function of that parameter, which
can be attributed to the different sharpness of the probability distribution.
We finally conclude by collecting and discussing our main results.

## The Model

### Orbifold $S^{1}/\mathbb{Z}_{2}$ and Background Geometry

As described in the Introduction, we consider the Randall-Sundrum model (RSM)
[28] with a single extra dimension $y$, together with a 2-dimensional ordinary
spacetime, whose coordinates are denoted by $x^{\mu}=\\{t,x\\}$. The total
spacetime possesses $D=3$ dimensions. The extra dimension $y$ is compactified
on a circle of radius $R$, and subject to a $\mathbb{Z}_{2}$ symmetry. These
features are captured by the equivalences

$\displaystyle S^{1}:\;y$ $\displaystyle\sim y+2\pi R\leavevmode\nobreak\ ,$
(1) $\displaystyle\mathbb{Z}_{2}:\;y$ $\displaystyle\sim-y\leavevmode\nobreak\
,$ (2)

which define the orbifold $S^{1}/\mathbb{Z}_{2}$ describing this extra
dimension. Along the $y$ dimension, the orbifold is a finite segment with two
fixed points at $y=0$ and $y=\pi R\equiv L$. The RSM assumes that there is a
$(D-1)$-brane of ordinary dimensions at each fixed point, see Fig. 1 for a
sketch of the space configuration and the orbifold symmetries.

Figure 1: Schematic representation of the extra dimension in the Randall-
Sundrum model.

The matter fields are supposed to reside on the brane at $y=L$, which is
referred to as the “visible brane”, while the brane at $y=0$ is the “hidden
brane”. Both branes contribute to the bulk background geometry through their
tensions, or vacuum energies, $T_{\mathrm{vis}}$ and $T_{\mathrm{hid}}$
respectively [28, 29]. The total background action is

$S=\int_{-L}^{L}dy\int
dx^{\mu}\sqrt{|g|}\left(2\alpha\mathcal{R}-\Lambda\right)+S_{\mathrm{vis}}+S_{\mathrm{hid}}\leavevmode\nobreak\
,$ (3)

where the first term is the usual Einstein-Hilbert action of the total space,
with $\Lambda$ the bulk cosmological constant, $\alpha$ a constant and $|g|$
the absolute value of the metric determinant, while

$\displaystyle S_{\mathrm{vis}}$ $\displaystyle=-\int
dx^{\mu}\sqrt{|g_{\mathrm{vis}}|}T_{\mathrm{vis}}\leavevmode\nobreak\ ,$ (4)
$\displaystyle S_{\mathrm{hid}}$ $\displaystyle=-\int
dx^{\mu}\sqrt{|g_{\mathrm{hid}}|}T_{\mathrm{hid}}\leavevmode\nobreak\ ,$ (5)

are the action contributions of the branes tensions, with the induced metrics
$g_{\mathrm{vis}}(x^{\mu})=g(x^{\mu},y=L)$ and
$g_{\mathrm{hid}}(x^{\mu})=g(x^{\mu},y=0)$. To address the hierarchy problem,
the following metric was proposed

$ds^{2}=e^{-2A(y)}\eta_{\mu\nu}dx^{\mu}dx^{\nu}-dy^{2}\leavevmode\nobreak\ ,$
(6)

where $e^{-2A(y)}$ is a warp factor, a rapidly changing function along the
additional dimension, and $\eta_{\mu\nu}$ is the Minkowski metric with
signature $(+,-)$. The metric in Eq. (6) obeys Einstein’s equations that are
obtained from the action (3): We refer the reader to the Supplementary
Material for the standard computation particularized to this lower-dimensional
spacetime. We also show that, as a consequence of these equations, the
function in the exponent is given by

$A(y)=k|y|$ (7)

where $k$ is the so called _warp coefficient_.

### Fermions in the Randall-Sundrum model

We now focus on the study of spin $1/2$ fermions, whose evolution equation is
the Dirac equation in curved spacetime

$(i\gamma^{a}e_{a}^{\mu}D_{\mu}-m)\Psi=0\leavevmode\nobreak\ .$ (8)

The $\gamma^{a}$ are the Dirac gamma matrices in a local rest frame, and the
covariant derivative is

$D_{\mu}=\partial_{\mu}-\frac{i}{4}\omega_{\mu}^{ab}\sigma_{ab}\leavevmode\nobreak\
,\quad\text{with}\quad\sigma_{ab}=\frac{i}{2}[\gamma_{a},\gamma_{b}]\leavevmode\nobreak\
,$ (9)

where $\omega_{\mu}^{ab}$ is the spin connection. The _vierbeins_
$e_{a}^{\mu}$ allow to express the Dirac matrices in a rest frame, that is,
they perform a change of basis to a non-coordinate system in which the metric
becomes the Minkowski metric

$g_{\mu\nu}e_{a}^{\mu}e_{b}^{\nu}=\eta_{ab}\leavevmode\nobreak\ .$ (10)

Equation (8) defines the vector current

$j^{\mu}=\sqrt{|g|}e_{a}^{\mu}\overline{\Psi}\gamma_{a}\Psi\leavevmode\nobreak\
,$ (11)

whose conservation $\partial_{\mu}j^{\mu}=0$ imposes the normalization
condition

$\int dx^{\mu}\sqrt{|g|}e_{0}^{0}\Psi^{\dagger}\Psi=1\leavevmode\nobreak\ .$
(12)

In the case of 2 spatial dimensions, the Dirac equation (8) can be reduced,
after some algebra, to

$i\gamma^{a}\left[e_{a}^{\mu}\partial_{\mu}\Psi+\frac{1}{2\sqrt{|g|}}\partial_{\mu}(e_{a}^{\mu}\sqrt{|g|})\Psi\right]-m\Psi=0\leavevmode\nobreak\
,$ (13)

where the $\gamma_{a}$ matrices become Pauli matrices. A simple choice of the
vierbeins obeying relation (10) is

$e_{0}=(e^{A(y)},0,0)\quad e_{1}=(0,e^{A(y)},0)\quad
e_{2}=(0,0,1)\leavevmode\nobreak\ ,$ (14)

which yields the following expression for the Dirac equation

$i\partial_{t}\Psi=-i\gamma^{0}\gamma^{1}\partial_{x}\Psi-i\gamma^{0}\gamma^{2}\partial_{y}(e^{-A(y)}\Psi)+\gamma_{0}e^{-A(y)}m\Psi\leavevmode\nobreak\
.$ (15)

This expression can be rewritten in Hamiltonian form as

$i\partial_{t}\chi=\mathcal{H}\chi,$ (16)

with

$\mathcal{H}=-\frac{i}{2}\\{B^{x},\partial_{x}\\}-\frac{i}{2}\\{B^{y},\partial_{y}\\}+\gamma_{0}e^{-A(y)}m\leavevmode\nobreak\
,$ (17)

where the change of variable $\chi=e^{-A(y)/2}\Psi$ was performed, and we
defined

$B^{x}=\gamma^{0}\gamma^{1}\leavevmode\nobreak\ ,\quad
B^{y}=e^{-A(y)}\gamma^{0}\gamma^{2}\leavevmode\nobreak\ .$ (18)

The symbol $\\{\cdot,\cdot\\}$ represents the anticommutator of two operators.
There is some freedom in the choice of the gamma matrices. For convenience, we
choose

$\gamma^{0}=\sigma_{x}\leavevmode\nobreak\
,\quad\gamma^{1}=i\sigma_{y}\leavevmode\nobreak\
,\quad\gamma^{2}=i\sigma_{z}\leavevmode\nobreak\ .$ (19)

### Boundary Conditions for Fermionic Fields

The periodic condition (1) simply implies that the fermionic fields need also
to be periodic

$\chi(x^{\mu},y+2L)=\chi(x^{\mu},y)\leavevmode\nobreak\ ,$ (20)

but the $\mathbb{Z}_{2}$ needs a deeper consideration, since it has to leave
the fermionic action invariant. We can write the fermionic action as

$S_{F}=\int dx^{\mu}\int
dy\overline{\chi}(x^{\mu},y)\left(i\gamma^{\mu}\partial_{\mu}+i\gamma^{2}\partial_{y}e^{-A(y)}-e^{-A(y)m}\right)\chi(x^{\mu},y)\leavevmode\nobreak\
.$ (21)

which is extremized by the Dirac equation (15). Under the action of
$\mathbb{Z}_{2}$ it becomes

$S_{F}=\int dx^{\mu}\int
dy\overline{\chi}(x^{\mu},-y)\left(i\gamma^{\mu}\partial_{\mu}-i\gamma^{2}\partial_{y}e^{-A(-y)}-e^{-A(-y)}m\right)\chi(x^{\mu},-y)\leavevmode\nobreak\
.$ (22)

We have to find an operator $M$, defined as
$\chi(x^{\mu},-y)=M\chi(x^{\mu},y)$, so as to allow the action to remain
invariant. Action (22) then becomes

$S_{F}=\int dx^{\mu}\int
dy\overline{\chi}(x^{\mu},y)\gamma^{0}M\gamma^{0}\left(i\gamma^{\mu}\partial_{\mu}-i\gamma^{2}\partial_{y}e^{-A(y)}-e^{-A(y)}m\right)M^{\dagger}\chi(x^{\mu},y)\leavevmode\nobreak\
,$ (23)

and establishes the following restrictions for $M$ to keep the action (21)
invariant,

$\displaystyle\gamma^{0}M\gamma^{0}\gamma^{\mu}M^{\dagger}=\gamma^{\mu}\leavevmode\nobreak\
,$ (24)
$\displaystyle\gamma^{0}M\gamma^{0}\gamma^{2}M^{\dagger}=-\gamma^{2}\leavevmode\nobreak\
,$ (25)
$\displaystyle\gamma^{0}M\gamma^{0}M^{\dagger}=\mathbb{I}\leavevmode\nobreak\
,$ (26)

where the first 2 conditions come from the kinetic terms of the action, and
the last one arises from the mass term. There does not exist a solution for
$M$ that solves all conditions simultaneously, although $M=\eta\sigma_{z}$ is
a solution for the first 2, with $\eta=\pm 1$. This means that a constant mass
term is forbidden. In the following we restrict ourselves to the case where
the “bulk mass” $m$ vanishes. The action of the fermionic field is therefore

$\begin{split}S_{F}=&\int
dx^{\mu}dy\overline{\chi}(x^{\mu},y)\left(i\gamma^{\mu}\partial_{\mu}+i\gamma^{2}\partial_{y}e^{-A(y)}\right)\chi(x^{\mu},y)\leavevmode\nobreak\
,\end{split}$ (27)

and the fermionic field has to obey the boundary condition

$\chi(x^{\mu},-y)=\eta\sigma_{z}\chi(x^{\mu},y)\leavevmode\nobreak\ ,$ (28)

with $\eta=\pm 1$.

Figure 2: Plots of the probability distribution for the first four stationary
states, with positive energy and a value of $q=10$, for $kL=3$ on the left and
for $kL=7$ on the right.

### Stationary Solutions

In this model, the Dirac field satisfies a complicated equation, Eq. (15),
which is difficult to address even numerically. In order to obtain some
insight, we first look for stationary solutions, which are defined as the
eigenstates of the Hamiltonian. For $m=0$, and with our choice of the gamma
matrices, the Hamiltonian takes the form

$\mathcal{H}=-\sigma_{z}\hat{p}_{x}+\frac{\sigma_{y}}{2}\left(e^{-A(y)}\hat{p}_{y}+\hat{p}_{y}e^{-A(y)}\right)\leavevmode\nobreak\
,$ (29)

where $p_{k}=-i\partial_{k}$ is the momentum operator along the $k$ direction
($k=x,y$). The stationary states $\phi_{n}(x,y)$ corresponding to energy
$E_{n}$ satisfy

$\mathcal{H}\phi_{n}(x,y)=E_{n}\phi_{n}(x,y)\leavevmode\nobreak\ .$ (30)

It is convenient to introduce a Fourier transform on the ordinary dimension
$x$:

$\tilde{\phi}_{n}(q,y)=\int dxe^{-iqx}\phi_{n}(x,y)\leavevmode\nobreak\ ,$
(31)

since the field is free to move along this direction. We found the energies

$E_{n}=\pm\sqrt{q^{2}+\left(k\alpha_{n}\right)^{2}}\leavevmode\nobreak\ ,$
(32)

where

$\alpha_{n}=\frac{n\pi}{e^{kL}-1}\leavevmode\nobreak\ ,\leavevmode\nobreak\
n=0,1,\;\dots$ (33)

The eigenfunctions associated with this spectrum that satisfy the boundary
condition (28), for the particular case with $\eta=1$, are

$\displaystyle\tilde{\phi}_{n}^{\uparrow}(q,y)$
$\displaystyle=\sqrt{\frac{2k}{e^{kL}-1}}\frac{E_{n}+q}{\sqrt{(E_{n}+q)^{2}+(k\alpha_{n})^{2}}}e^{\frac{k|y|}{2}}\cos\left[\alpha_{n}\left(e^{k|y|}-1\right)\right]\leavevmode\nobreak\
,$ (34) $\displaystyle\tilde{\phi}_{n}^{\downarrow}(q,y)$
$\displaystyle=\sqrt{\frac{2k}{e^{kL}-1}}\frac{k\alpha_{n}}{\sqrt{(E_{n}+q)^{2}+(k\alpha_{n})^{2}}}e^{\frac{k|y|}{2}}\sin\left[\alpha_{n}\left(e^{k|y|}-1\right)\right]\mathrm{sign}(y)\leavevmode\nobreak\
,$ (35)

where the components of the spinor field are
$\tilde{\phi}_{n}=(\tilde{\phi}_{n}^{\uparrow},\tilde{\phi}_{n}^{\downarrow})^{T}$.
The particular case $n=0$ only has an upper component, which is given by

$\tilde{\phi}_{0}^{\uparrow}(q,y)=\sqrt{\frac{k}{e^{kL}-1}}e^{\frac{k|y|}{2}}\mathrm{sign}(E_{n}+q)\leavevmode\nobreak\
,$ (36)

and is undefined for energy and momentum with different sign. The procedure to
obtain the eigenfunctions is detailed in the Supplementary Material, as well
as the solution for $\eta=-1$. The probability distribution associated to
these wavefunctions is concentrated around $y=L$ for high values of the warp
coefficient $k$. We illustrate this behavior in Fig. 2, where we have plotted
the probability density for the first modes with positive energy, and momentum
$q=10$, for a value of the warp coefficient $kL=3$ and $kL=7$, respectively.

## A Quantum Walk for the Randall-Sundrum model

Once we have discussed the main properties of the RSM in the continuum
spacetime, we focus on the main goal of our work, which consists in
constructing a QW that is able to simulate the dynamics of a spin $1/2$
particle subject to the geometric effects and symmetries of the model. To
incorporate the metric, we adapt the scheme introduced in [9], which allows to
reproduce (in the continuum limit) a Dirac equation of the form Eq. (16). The
QW is defined on a 2-dimensional discrete grid with $x$ and $y$ axis, with
discrete positions labeled by $r$ and $s$, respectively. The grid points are
equally spaced by $\epsilon$, so that the spatial coordinates can be related
to the grid points by $x=\epsilon r$ and $y=\epsilon s$. The Hilbert space
that corresponds to these spatial degrees of freedom, $\mathcal{H}_{spatial}$
is spanned by the basis $\\{|x=\epsilon r,y=\epsilon
s\rangle\\}/r,s\in\mathbb{Z}$. Time steps are labeled by $j\in\mathbb{Z}$, and
are also equally spaced by $\epsilon$. The coin (or internal) space is a 2
dimensional Hilbert space $\mathcal{H}_{\mathrm{coin}}$, so that the total
Hilbert space is
$\mathcal{H}_{tot}=\mathcal{H}_{spatial}\otimes\mathcal{H}_{\mathrm{coin}}$.
At a given time step, the state of the walker will be represented by a two
component spinor $|\chi_{j}\rangle\in\mathcal{H}_{tot}$. The one step
evolution of the QW is given by

$|\chi_{j+1}\rangle=U|\chi_{j}\rangle\leavevmode\nobreak\ ,$ (37)

where $U$ is a unitary operator which consists on rotations in the components
of $|\chi_{j}\rangle$, and translations by $\epsilon$ in the two directions of
physical space $x$ and $y$. The angles of rotation can, in general, be
dependent on the spacetime coordinates of the walker. Following [9], we adopt

$U=R^{-1}(y)\left[\Theta(y)S_{y}(\epsilon/2)\right]^{2}R(y)S_{x}(\epsilon)\leavevmode\nobreak\
,$ (38)

where $S_{k}(\epsilon)=\exp(-i\sigma_{z}p_{k}\epsilon)$ are spin-dependent
shift operators in the direction $\pm k$ (with $k=x,y$),

$\Theta(y)=\begin{pmatrix}c(y)&is(y)\\\
-is(y)&c(y)\end{pmatrix}\leavevmode\nobreak\ ,$ (39)

with $c(y)=e^{-A(y)}$, $s(y)=\sqrt{1-e^{-2A(y)}}$, and

$R(y)=\frac{1}{\sqrt{2}}\begin{pmatrix}f^{*}(y)&if(y)\\\
-f^{*}(y)&-if(y)\end{pmatrix}\leavevmode\nobreak\ ,$ (40)

where $f(y)=\sqrt{\frac{1+c(y)}{2}}+i\sqrt{\frac{1-c(y)}{2}}$. At each
position $(r,s)$ we introduce

$\chi_{j,r,s}\equiv\langle x=\epsilon r,y=\epsilon
s|\chi_{j}\rangle=\begin{pmatrix}\chi_{j,r,s}^{\uparrow}\\\
\chi_{j,r,s}^{\downarrow}\end{pmatrix}\leavevmode\nobreak\ ,$ (41)

which represents the amplitude (given a component of the spin) for the
particle to be localized at the position labeled by $(r,s)$ and time step $j$.
In this way, the time step defined by (37) can be recast as a recursive
formula for $\chi_{j,r,s}$, which is provided in the Supplementary Material.
In order to implement this QW to simulate fermions in the RSM, appropriate
conditions have to be set to comply with the boundary conditions (20) and
(28). It can be explicitly shown, from the recursive formula for
$\chi_{j,r,s}$, that this QW dynamics respects (28), in the sense that, if the
walker obeys the condition

$\chi_{j,r,-s}=\eta\sigma_{z}\chi_{j,r,s}\leavevmode\nobreak\ ,$ (42)

at time $j$, it is also obeyed at time $j+1$. For the simulations, we
discretize the $y$ coordinate along the segment $[-L,L]$ with a spacing
$\epsilon$, and impose an initial condition which satisfies Eq. (42). We use
the same lattice spacing in the $x$ direction, together with an strategy that
adapts its effective extension to the time step. We also impose periodic
boundary conditions on the grid to respect condition (20), taking into account
that functions evaluated at $y=L+\epsilon$ should be identified with functions
at $y=-L+\epsilon$ to respect the periodicity in the range $[-L,L]$.

## Results

The QW defined in the previous section is guaranteed to reproduce (in the
continuum limit) a Dirac equation of the form (16), such as the one
corresponding to the RSM. The question that arises concerns the dynamics
appearing at a finite lattice and time step spacing. Of course, one does not
expect the QW to behave exactly as the continuum field but, to what extent do
they differ? Are there any new features that appear in the discrete case? In
particular, we are interested in looking for some kind of probability
concentration towards the visible brane, for a given initial condition. In
this Section we explore all these features.

### Stationarity of the Eigenstates Solutions on the Quantum Walk

As an initial comparison, we start by considering the discretized version of
the eigenstates corresponding to the continuum limit Hamiltonian, obtained
before. Such states remain stationary within this limit (i.e. they just evolve
by adopting a trivial phase). How do they evolve under the action of the QW?
We consider an initial state which corresponds to an eigenstate of the
continuum, with fixed momentum $q$, and check whether the QW evolution of this
state is stationary. The initial condition of the walker is therefore

$\chi_{0,r,s}=\tilde{\phi}_{n}(q,\epsilon s)e^{iq\epsilon
r}\leavevmode\nobreak\ ,$ (43)

which represents a constant probability density along the ordinary dimension
$x$. As expected, the QW evolution does not remain stationary, although it
keeps a close resemblance to the initial state. This can be observed from Fig.
3, where we represented the normalized marginal probability along the $y$
direction of the walker (after summing over $x$) at different time steps, for
an initial stationary state solution with $n=2$, and warp coefficient $kL=3$.

Figure 3: Snapshots of the probability density starting from an initial
eigenstate with $n=2$ and positive energy, for a value of $kL=3$, and $q=10$.
The simulation grid has 100 points along the $y$ direction, and enough points
have been taken in the $x$ direction to ensure that the total probability
density does not leak outside the boundaries.

### Localization in the QW

We now investigate the localization capability of the above defined QW, i.e.,
whether it shows a tendency to concentrate the walker towards the visible
brane at $y=L$. We consider an initial walker which is fully localized

$\chi_{0,r,s}=\delta_{x,0}\delta_{y,y_{0}}C_{0}\leavevmode\nobreak\ ,$ (44)

where $C_{0}$ is the initial coin state, and we recall that $x=\epsilon r$ and
$y=\epsilon s$. We explore the evolution of a walker which is initially
localized at the center of the extra dimension, that is at
$y_{0}=\frac{L}{2}$, and we study the probability distribution for different
values of the warp coefficient, at a given time step. In Fig. 4 we show the
surface plot of the probability density with the above initial conditions, and
$C_{0}=\frac{1}{\sqrt{2}}(1,i)^{T}$, which induces a symmetric evolution in
the ordinary dimension. The blue (red) color of the surface represents
dominance of the upper (lower) coin component, while yellow stands for a
superposition of both components.

Figure 4: Probability density distribution, at $t=5L$, of an initial localized
walker centred at $(x_{0},y_{0})=(0,L/2)$ for different values of $kL$ with
initial coin components $C_{0}=\frac{1}{\sqrt{2}}(1,i)^{T}$. The height of the
curve represents the probability of finding the walker in that position, and
the colors indicate the coin state. The red (blue) color indicates a
predominance of the upper (lower) component, while yellow stands for a
superposition of both components. The simulation grid has 100 points along the
$y$ direction, and enough points have been taken in the $x$ direction to
ensure that the total probability density does not leak outside the
boundaries. Figure 5: Expected value of the probability distribution along the
extra dimension $y$, as calculated from the FPD, for different values of the
warp coefficient $kL$. The initial condition is the same as in Fig. 4.

We notice that most of the probability distribution in the $x$ direction is
concentrated along a freely propagating front which moves at the maximum speed
($x=\pm t$), consistently with the fact that the QW simulates massless
fermions. We also notice that most of the right propagating distribution
(positive values of $x$) is dominated by the upper coin component, while the
part propagating to the left (negative values of $x$) mainly contains the
lower coin component, a fact that can also be inferred from the explicit
evolution of the QW (see Supplementary Material for details). The propagation
of the walker along the extra dimension $y$ strongly depends on the value of
the warp coefficient. At $t=5L$, the distribution with the lowest value of
$kL$ possesses non-zero values on the visible brane $y=L$, while the other two
do not. In fact, the displacement of the probability distribution towards
$y=L$ is slower for the highest $kL$. In other words, a larger value of the
warp coefficient dramatically increases the time scale of the dynamics along
the extra dimension, and makes it prohibitively expensive (in terms of
computational cost) to explore larger values of $kL$ than those considered
here.

In order to investigate whether the QW exhibits the same behavior as the
stationary states, in the sense that a higher value of the warp coefficient
induces a stronger localization near the visible brane, we study the
distribution of the freely propagated parts of the walker (the regions around
$x=\pm t$), where most of the probability density is concentrated, as can be
readily seen in Fig. 4. The probability distribution associated to these two
zones will be referred to as the “freely propagating distribution” (FPD). In
terms of the spinor components, those are the probability density
distributions obtained from $\chi_{j,s}^{\mathrm{R}}\equiv\chi_{j,j,s}$ and
$\chi_{j,s}^{\mathrm{L}}\equiv\chi_{j,-j,s}$, where $r=\pm j$ restricts the
wavefunction to the two freely propagating peaks. In Fig. 5 we represent the
expected value for these distributions along the $y$ dimension, which can be
defined as

$\langle y_{R(L)}(t)\rangle=\sum_{s}\epsilon
s\;\chi_{j,s}^{\mathrm{R(L)}\dagger}\chi_{j,s}^{\mathrm{R(L)}}\leavevmode\nobreak\
,$ (45)

where $t=\epsilon j$, for different values of $kL$. First of all we notice
that this quantity reaches an asymptotic value, which is closer to $L$ for
higher warp coefficients. Secondly, as discussed above, the warp coefficient
induces a change in the time scale of the dynamics, so that lower values of
the warp coefficient show a faster convergence towards the asymptotic state,
consistently with the features already observed in Fig. 4.

Figure 6: Probability distributions of the FPDs along the extra dimension $y$,
for the value $kL=3$. The inset is a histogram showing the value of the
$B_{n}(t)$ coefficients, as defined by Eq. (47): see the text for an
explanation. The left (right) panels show the left (right) FPD. The top panels
are calculated at a shorter time $t=50L$ and the bottom ones at a longer time
$t=1000L$. The initial condition is the same as in Fig. 5, and the simulation
grid has 200 points along the $y$ direction.

### Mode decomposition of the freely propagating distribution

Our simulations indicate that the FPD reaches a steady state along the extra
dimension, in a similar fashion as the expected value (45). This evolution can
be appreciated from the plots of Fig. 6. Al late times (lower row), the
probability distribution resembles the probability density of a stationary
state with positive energy and momentum in one of the lowest modes: $n=0$ for
the right FPD, and $n=1$ for the left FPD. It is important to recall that, as
discussed above, the right (left) FPD is predominantly composed by the upper
(lower) component of the spinor, and that $n=0$ has no lower component: see
Eq. (36). This causes a fundamental difference when comparing the left and
right contributions. In order to investigate these features on the time
evolution, we introduce a decomposition on the wavefunction of the walker as a
combination of the stationary states basis. This allows us to write

$\chi_{j,r,s}=\int_{-\pi/\epsilon}^{\pi/\epsilon}\frac{dq}{2\pi}\sum_{n}\beta_{n}(q,t)\tilde{\phi}_{n}(q,\epsilon
s)e^{-iq\epsilon r}\leavevmode\nobreak\ ,$ (46)

where the temporal dependence is included on the $\beta_{n}(q,t)$
coefficients. In the Supplementary Material we detail how these factors can be
computed, and define their normalization conditions. In particular, we are
interested on the contribution of each value $n$, therefore we integrate out
the dependence in the quasi-momentum $q$. In other words, we are interested on
the following (time-dependent) coefficients:

$B_{n}(t)=\int_{-\pi/\epsilon}^{\pi/\epsilon}\frac{dq}{2\pi}\left|\beta_{n}(q,t)\right|^{2}\leavevmode\nobreak\
.$ (47)

The different mode components $B_{n}(t)$ of Fig. 6 have been included as an
inset in those plots. On the one hand, it can be observed that, at long times,
when a steady state has been reached, the FPDs are mostly composed by the
lowest possible mode ($n=0$ or $n=1$, as discussed above). On the other hand,
at short times, the FPDs contain additional higher modes.

### Entanglement Entropy

Finally, we study the entanglement properties that the QW exhibits between the
coin and position degrees of freedom for the already considered, initially
localized state. The entanglement can be quantified using the von Neumann
entropy of the reduced density matrix in the coin space

$S(t)=-\mathrm{Tr}\left\\{\rho_{c}(t)\log_{2}\rho_{c}(t)\right\\}\leavevmode\nobreak\
,$ (48)

where $\rho_{c}(t=\epsilon j)=\sum_{r,s}\chi_{j,r,s}\chi_{j,r,s}^{\dagger}$ is
the reduced density matrix in the coin space, i.e. after tracing out the
spatial degrees of freedom. In Fig. 7 we plot the evolution of the
entanglement entropy of a fully localized initial state for different values
of the warp coefficient, with a coin state
$C_{0}=\frac{1}{\sqrt{5}}(1,2i)^{T}$. Notice that this choice is different
from that one used in the previous section, for reasons that are explained
below. It can be seen that the entanglement entropy reaches lower values as
$kL$ increases, an effect that can probably be due to the fact that the
probability density in between the FPDs becomes more spread (and therefore
“less ordered”) at lower values of $kL$. This can be observed in Fig. 8, where
we plotted a zoomed version of Fig. 4, but obtained with the above initial
coin components $C_{0}=\frac{1}{\sqrt{5}}(1,2i)^{T}$. One can see that, for
lower values of the warp coefficient, a significant part of the probability
distribution is scrambled in the intermediate region between both parts of the
FPD. This diffusion effect can be totally mitigated for extreme values of the
warp coefficient, leading to a minimum value of the entropy which is
completely dominated by the FPD, and can be obtained from the initial coin
components. In the Supplementary Material we show this limiting situation, and
how the corresponding entropy can be computed. The initial coin state
$C_{0}=\frac{1}{\sqrt{2}}(1,i)^{T}$ previously used produces values of the
entropy which are very close to unity in all cases, making it difficult to
appreciate the effects that are discussed above.

Figure 7: Evolution of the entanglement entropy with the initial condition Eq.
(44) centred at $y_{0}=L/2$ for different values of the warp coefficient and
initial coin components $C_{0}=\frac{1}{\sqrt{5}}(1,2i)$. The dotted line
represent the minimum value the entropy can reach for very high values of
$kL$, which is computed in the Supplementary Material. This simulation grid
has 50 points along the $y$ direction. Figure 8: Probability density
distribution, at $t=5L$, for an initially localized walker centered at
$(x_{0},y_{0})=(0,L/2)$, and different values of $kL$, with initial coin
components $C_{0}=\frac{1}{\sqrt{5}}(1,2i)^{T}$. The vertical axis has been
zoomed in to show that the probability density between the two regions of the
FPD is more scrambled for lower values of the warp coefficient. The colors and
grid parameters are the same as in Fig. 4.

## Conclusions

We have investigated a quantum walk which allows to simulate the Randall-
Sundrum model of extra dimensions, while satisfying the constrains imposed by
the symmetries of that model. This model has played an important role in high
energy physics, aiming to solve the hierarchy problem, by introducing one
finite extra dimension that possesses two branes at its extremes. The matter
fields are confined in the visible brane, while gravity is allowed to span
along this whole dimension. We worked it out for the case of spin $1/2$
fermions in a two dimensional space, composed by an ordinary dimension and an
orbifolded one, apart from time, and obtained the Dirac equation in this
spacetime configuration. The boundary conditions of the orbifold on the
fermionic field forced it to be massless on the bulk. In this lower
dimensional space we were able to obtain the eigenenergies of the fermionic
field, as well as the corresponding eigenstates, showing a probability density
which is concentrated near the visible brane, a phenomenon that bears an
analogy with the localization effect that can be found in many scenarios [19,
20, 21, 23, 24, 26, 27].

This analogy motivated us to seek localization effects on the QW that we
introduced to simulate the RSM. The QW is defined in such a way that, in the
continuum limit, the Dirac equation of the fermionic field for the RSM metric
is recovered. We investigated the confining capabilities of the QW, by
considering an initially localized walker away from the visible brane. We
concluded that the freely propagating parts of the probability distribution,
where the probability is mostly concentrated, reach an asymptotic value of the
expected position along the extra dimension. Moreover, the asymptotic value
gets closer to the visible brane for higher values of the warp coefficient,
which therefore drives the strength of localization, and also noticed that it
had an effect on the timescale of the dynamics, by delaying them for higher
values of the coefficient.

At long time steps, the probability densities show an asymptotic shape, with a
resemblance with the eigenstates that were obtained in the continuous model,
which suggested a study based on the decomposition of the wavefunction in
terms of these stationary states. We found that the freely propagating parts
of the QW are dominated, in the asymptotic regime, by the lowest possible
(i.e., compatible with the symmetries of the model) modes. At intermediate
time steps, the same decomposition manifests a combination of multiple modes
with higher energy.

Finally, we found that the entanglement between coin and spatial degrees of
freedom is reduced for stronger warp coefficients. We associated this result
to the higher spreading of the density distribution for the lower values of
the warp coefficient.

We conclude that quantum walks are suitable candidates for simulating models
of field theories with extra dimensions that rely on the curvature of the
spacetime. Not only the model is interesting from the point of view of the
field theory: It allows to design a quantum process that can be tailored to
exhibit very rich dynamics, showing free propagation in one dimension, and an
asymptotic confining behavior on the other one, with rates that can be tuned
by an appropriate choice of the parameters. In this way, the interplay between
high energy physics and quantum simulations can be of mutual benefit.

## Acknowledgements

This work has been founded by the Spanish FEDER/MCIyU-AEI grant
FPA2017-84543-P, SEV-2014-0398 and Generalitat Valenciana grant
PROMETEO/2019/087. We also acknowledge support from CSIC Research Platform
PTI-001.

## References

  * [1] Cedzich, C. _et al._ Propagation of quantum walks in electric fields.  _Phys. Rev. Lett._ 111, 160601, DOI: 10.1103/PhysRevLett.111.160601 (2013).
  * [2] Arnault, P. & Debbasch, F.  Landau levels for discrete-time quantum walks in artificial magnetic fields.  _Physica A: Statistical Mechanics and its Applications_ 443, 179 – 191, DOI: https://doi.org/10.1016/j.physa.2015.08.011 (2016).
  * [3] Di Molfetta, G., Brachet, M. & Debbasch, F.  Quantum walks in artificial electric and gravitational fields.  _Physica A: Statistical Mechanics and its Applications_ 397, 157–168 (2014).
  * [4] Yalç ınkaya, i. d. I. & Gedik, Z.  Two-dimensional quantum walk under artificial magnetic field.  _Phys. Rev. A_ 92, 042324, DOI: 10.1103/PhysRevA.92.042324 (2015).
  * [5] Bru, L. A., Hinarejos, M., Silva, F., de Valcárcel, G. J. & Roldán, E.  Electric quantum walks in two dimensions.  _Phys. Rev. A_ 93, 032333, DOI: 10.1103/PhysRevA.93.032333 (2016).
  * [6] Arnault, P. & Debbasch, F.  Quantum walks and discrete gauge theories.  _Phys. Rev. A_ 93, 052301, DOI: 10.1103/PhysRevA.93.052301 (2016).
  * [7] Márquez-Martín, I., Arnault, P., Di Molfetta, G. & Pérez, A.  Electromagnetic lattice gauge invariance in two-dimensional discrete-time quantum walks.  _Phys. Rev. A_ 98, 032333, DOI: 10.1103/PhysRevA.98.032333 (2018).
  * [8] Cedzich, C., Geib, T., Werner, A. H. & Werner, R. F.  Quantum walks in external gauge fields.  _Journal of Mathematical Physics, Vol.60, Issue 1, 2019_ DOI: 10.1063/1.5054894 (2019).  1808.10850v1.
  * [9] Arnault, P. & Debbasch, F.  Quantum walks and gravitational waves.  _Annals of Physics_ 383, 645 – 661, DOI: http://dx.doi.org/10.1016/j.aop.2017.04.003 (2017).
  * [10] Arrighi, P., Di Molfetta, G., Marquez-Martin, I. & Perez, A.  From curved spacetime to spacetime-dependent local unitaries over the honeycomb and triangular quantum walks.  _Scientific Reports_ 9, 10904 (2019).
  * [11] Molfetta, G. D. & Pérez, A.  Quantum walks as simulators of neutrino oscillations in a vacuum and matter.  _New Journal of Physics_ 18, 103038, DOI: 10.1088/1367-2630/18/10/103038 (2016).
  * [12] Mallick, A., Mandal, S. & Chandrashekar, C. M.  Neutrino oscillations in discrete-time quantum walk framework.  _The European Physical Journal C_ 77, 85 (2017).
  * [13] Jha, A. K., Chatla, A. & Bambah, B. A.  Quantum simulation of oscillating neutrinos.  _5th International Conference on Particle Physics and Astrophysics_ (2020).  2010.06458v2.
  * [14] Arnault, P., Perez, A., Arrighi, P. & Farrelly, T.  Discrete-time quantum walks as fermions of lattice gauge theory.  _Phys. Rev. A_ 99, 032110, DOI: 10.1103/PhysRevA.99.032110 (2019).
  * [15] Kaluza, T.  Zum unitätsproblem der physik.  _Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.)_ K1, 966, DOI: 10.1142/S0218271818700017 (1921).
  * [16] Klein, O.  Quantentheorie und fünfdimensionale relativitätstheorie.  _Zeitschrift für Physik_ 37, 895–906, DOI: 10.1007/BF01397481 (1926).
  * [17] Rubakov, V. & Shaposhnikov, M.  Do we live inside a domain wall?  _Physics Letters B_ 125, 136 – 138, DOI: http://dx.doi.org/10.1016/0370-2693(83)91253-4 (1983).
  * [18] Márquez-Martín, I., Di Molfetta, G. & Pérez, A.  Fermion confinement via quantum walks in (2+1)-dimensional and (3+1)-dimensional space-time.  _Phys. Rev. A_ 95, 042112, DOI: 10.1103/PhysRevA.95.042112 (2017).
  * [19] Anderson, P. W.  Absence of Diffusion in Certain Random Lattices.  _Phys. Rev._ 109 (1956).
  * [20] Aubry, S. & André, G.  Analyticity breaking and anderson localization in incommensurate lattices.  _Ann. Israel Phys. Soc_ 3, 18 (1980).
  * [21] Grempel, D. R., Fishman, S. & Prange, R. E.  Localization in an incommensurate potential: An exactly solvable model.  _Phys. Rev. Lett._ 49, 833–836, DOI: 10.1103/PhysRevLett.49.833 (1982).
  * [22] Lahini, Y. _et al._ Observation of a localization transition in quasiperiodic photonic lattices.  _Phys. Rev. Lett._ 103, 013901, DOI: 10.1103/PhysRevLett.103.013901 (2009).
  * [23] Joye, A. & Merkli, M.  Dynamical localization of quantum walks in random environments.  _Journal of Statistical Physics_ 140, 1025–1053, DOI: 10.1007/s10955-010-0047-0 (2010).
  * [24] Schreiber, A. _et al._ Decoherence and disorder in quantum walks: From ballistic spread to localization.  _Phys. Rev. Lett._ 106, 180403, DOI: 10.1103/PhysRevLett.106.180403 (2011).
  * [25] Crespi, A. _et al._ Anderson localization of entangled photons in an integrated quantum walk.  _Nature Photonics_ 7, 322–328, DOI: 10.1038/nphoton.2013.26 (2013).  arXiv:1304.1012v1.
  * [26] Navarrete-Benlloch, C., Pérez, A. & Roldán, E.  Nonlinear optical Galton board.  _Physical Review A - Atomic, Molecular, and Optical Physics_ 75, 1–7, DOI: 10.1103/PhysRevA.75.062333 (2007).  0604084.
  * [27] Shikano, Y. & Katsura, H.  Localization and fractality in inhomogeneous quantum walks with self-duality.  _Phys. Rev. E_ 82, 031122, DOI: 10.1103/PhysRevE.82.031122 (2010).
  * [28] Randall, L. & Sundrum, R.  Large mass hierarchy from a small extra dimension.  _Phys. Rev. Lett._ 83, 3370–3373, DOI: 10.1103/PhysRevLett.83.3370 (1999).
  * [29] Arkani-Hamed, N., Dimopoulos, S. & Dvali, G.  Phenomenology, astrophysics, and cosmology of theories with submillimeter dimensions and tev scale quantum gravity.  _Phys. Rev. D_ 59, 086004, DOI: 10.1103/PhysRevD.59.086004 (1999).

## Supplementary Material

### Metric Solution

Here we prove that the metric (6) extremizes the background action (3), which
can be rewritten as

$S=\int dy\int
dx^{\mu}\sqrt{|g|}\Big{(}2\alpha\mathcal{R}-\Lambda-\delta(y)T_{\mathrm{hid}}-\delta(y-L)T_{\mathrm{vis}}\Big{)}\leavevmode\nobreak\
.$ (49)

The extrema of this action gives the following Einstein equations

$\sqrt{|g|}\bigg{(}R_{MN}-\frac{1}{2}g_{MN}\mathcal{R}+\frac{1}{4\alpha}\Lambda
g_{MN}\bigg{)}=-\frac{\sqrt{|g|}}{4\alpha}\bigg{(}T_{\mathrm{hid}}\delta(y)g_{\mu\nu}\delta_{M}^{\mu}\delta_{N}^{\nu}+T_{\mathrm{hid}}\delta(y-L)g_{\mu\nu}\delta_{M}^{\mu}\delta_{N}^{\nu}\bigg{)}\leavevmode\nobreak\
,$ (50)

with indices $M,N=\\{t,x,y\\}$, while $\mu,\nu=\\{t,x\\}$ only account for
ordinary dimensions. After computing the curvature tensor, we obtain the
following equation for the $yy$ component

$A^{\prime}(y)^{2}+\frac{\Lambda}{4\alpha}=0\leavevmode\nobreak\ ,$ (51)

which yields the solution

$A(y)=k|y|\leavevmode\nobreak\ ,\quad\mathrm{with}\quad
k\equiv\sqrt{-\frac{\Lambda}{4\alpha}}\leavevmode\nobreak\ ,$ (52)

and has to be consistent with the orbifold symmetry (2). There are no $\mu y$
components, as the metric and tensors with these components vanish. The
$\mu\nu$ components of the Einstein equations are

$(A^{\prime}(y)^{2}-A^{\prime\prime}(y))e^{-2A(y)}\eta_{\mu\nu}+\frac{1}{4\alpha}\Lambda
e^{-2A(y)}\eta_{\mu\nu}=-\frac{1}{4\alpha}e^{-2A(y)}\eta_{\mu\nu}\left[T_{\mathrm{hid}}\delta(y)+T_{\mathrm{vis}}\delta(y-L)\right]\leavevmode\nobreak\
,$ (53)

which, making use of equation (51), can be simplified to

$A^{\prime\prime}(y)=\frac{1}{4\alpha}\left[T_{\mathrm{hid}}\delta(y)+T_{\mathrm{vis}}\delta(y-L)\right]\leavevmode\nobreak\
.$ (54)

After computing the second derivative of $A(y)$ from Eq. (52), and taking into
account the periodicity of the metric (1), yields

$A^{\prime\prime}(y)=2k(\delta(y)-\delta(y-L))\leavevmode\nobreak\ ,$ (55)

which allows us to identify, from Eq. (54), the values of the tensions

$-T_{\mathrm{vis}}=T_{\mathrm{hid}}=8\alpha
k=\sqrt{-16\Lambda\alpha}\leavevmode\nobreak\ .$ (56)

The results obtained in this section indicate that the bulk geometry has to be
Anti-de Sitter, with a negative bulk cosmological constant, and that the
visible brane has negative tension, while the hidden one is positive. These
results differ, from standard works on Randall-Sundrum, on the constant
coefficient appearing in the expression for $k$, Eq. (52) because we are
considering a one dimensional ordinary space, so that the computations of the
curvature tensor yield different constant factors.

### Hamiltonian eigenstates

In order to solve the eigenvalue problem (30), it is convenient to perform the
change of basis $\xi(q,y)=H\tilde{\phi}(q,y)$, with

$H=\begin{pmatrix}1&1\\\ 1&-1\end{pmatrix}$ (57)

the Hadamard matrix, so that the eigenvalue equation becomes

$\displaystyle\left[q+\frac{i}{2}\left(p_{y}e^{-A(y)}+e^{-A(y)}p_{y}\right)\right]\xi_{n}^{-}=E_{n}\xi_{n}^{+}\leavevmode\nobreak\
,$ (58)
$\displaystyle\left[q-\frac{i}{2}\left(p_{y}e^{-A(y)}+e^{-A(y)}p_{y}\right)\right]\xi_{n}^{+}=E_{n}\xi_{n}^{-}\leavevmode\nobreak\
,$ (59)

where $\xi^{\pm}$ are the components of $\xi=(\xi^{+},\xi^{-})^{T}$. This
system of equations can be decoupled, giving

$\left[q^{2}+\frac{1}{4}\left(p_{y}e^{-A(y)}+e^{-A(y)}p_{y}\right)^{2}\right]\xi_{n}^{\pm}=E_{n}^{2}\xi_{n}^{\pm}\leavevmode\nobreak\
,$ (60)

which is a second order differential equation that can be solved for the
appropriate boundary conditions. We solve this equation both in the positive
$\left[\xi^{\pm}(0<y<L)\right]_{P}$ and negative domain
$\left[\xi^{\pm}(-L<y<0)\right]_{N}$, delivering

$\displaystyle\left[\xi_{n}^{\pm}(y)\right]_{P}=$ $\displaystyle
Ae^{\frac{ky}{2}}\cos\left(e^{ky}\alpha_{n}\right)+Be^{\frac{ky}{2}}\sin\left(e^{ky}\alpha_{n}\right)\leavevmode\nobreak\
,$ (61) $\displaystyle\left[\xi_{n}^{\pm}(y)\right]_{N}=$ $\displaystyle
Ce^{-\frac{ky}{2}}\cos\left(e^{-ky}\alpha_{n}\right)+De^{\frac{ky}{2}}\sin\left(e^{-ky}\alpha_{n}\right)\leavevmode\nobreak\
,$ (62)

where we defined

$\alpha_{n}=\frac{\sqrt{E_{n}^{2}-q^{2}}}{k}\leavevmode\nobreak\ .$ (63)

These solutions are related by the continuity conditions

$\left[\xi_{n}^{\pm}(0)\right]_{P}=\left[\xi_{n}^{\pm}(0)\right]_{N}\leavevmode\nobreak\
,\;\text{and}\;\left[\xi_{n}^{\pm}(L)\right]_{P}=\left[\xi_{n}^{\pm}(-L)\right]_{N}\leavevmode\nobreak\
,$ (64)

where in the last one the periodicity of the wavefunctions, Eq. (20), has been
used, and imply that the solutions are related by $A=C$ and $B=D$. The
discontinuity introduced by the delta terms at $y=0$ and $y=\pm L$, coming
from $A^{\prime\prime}(y)$, imposes $B=A\tan\alpha_{n}$, and the following
restrictions to the energies

$\tan\alpha_{n}=\tan\left(e^{kL}\alpha_{n}\right)\leavevmode\nobreak\ ,$ (65)

which yields the spectrum in Eq. (32). After taking into account these
conditions, the eigenstates become

$\xi_{n}^{\pm}(y)=Ae^{\frac{k|y|}{2}}\left[\cos\left(e^{k|y|}\alpha_{n}\right)+\tan\alpha_{n}\sin\left(e^{k|y|}\alpha_{n}\right)\right]\leavevmode\nobreak\
.$ (66)

However, these solutions come from the second order differential equation
(60), whereas the original equations were first order, and relate $\xi^{+}(y)$
to $\xi^{-}(y)$. To find the appropriate solution of the eigenfunctions, we
need to take into account these relations. Since any lineal combination of
solutions is also a solution of the equations, we consider the solution
$\left[\xi_{n}^{+}(y)\right]_{2}$ of Eq. (60) to obtain
$\left[\xi_{n}^{-}(y)\right]_{1}$ from Eq. (59), where $[\leavevmode\nobreak\
\cdot\leavevmode\nobreak\ ]_{i}$ denotes whether the solution comes from a
first ($i=1$) or second ($i=2$) order differential equation. Similarly, from
$\left[\xi_{n}^{-}(y)\right]_{2}$ we obtain $\left[\xi_{n}^{+}(y)\right]_{1}$,
so that

$\left[\xi_{n}^{\pm}(y)\right]_{1}=\bigg{[}\sin\left(\alpha_{n}e^{k|y|}\right)\left(\frac{q}{E}\tan\alpha_{n}\pm\frac{k\alpha_{n}}{E}\mathrm{sign}(y)\right)+\cos\left(\alpha_{n}e^{k|y|}\right)\left(\frac{q}{E}\mp\frac{k\alpha_{n}}{E}\tan\alpha_{n}\mathrm{sign}(y)\right)\bigg{]}Ae^{\frac{k|y|}{2}}\leavevmode\nobreak\
.$ (67)

The general solution for the eigenstates is a lineal combination of this pair
of solutions

$\displaystyle\xi_{n}^{+}(y)=K_{1}\left[\xi_{n}^{+}(y)\right]_{2}+K_{2}\left[\xi_{n}^{+}(y)\right]_{1}\leavevmode\nobreak\
,$ (68)
$\displaystyle\xi_{n}^{-}(y)=K_{1}\left[\xi_{n}^{-}(y)\right]_{1}+K_{2}\left[\xi_{n}^{-}(y)\right]_{2}\leavevmode\nobreak\
,$ (69)

where the relation between the constants $K_{1}$ and $K_{2}$ is set by Eq.
(28), which, depending on the possible values of $\eta$, implies the
restrictions

$\displaystyle\eta=+1$ $\displaystyle\implies K_{1}=K_{2}\leavevmode\nobreak\
,$ (70) $\displaystyle\eta=-1$ $\displaystyle\implies
K_{1}=-K_{2}\leavevmode\nobreak\ .$ (71)

Finally, undoing the change of basis, we recover the original eigenstate
components of Eqs. (34,35) for $\eta=+1$, while

$\displaystyle\phi_{n}^{\uparrow}(y)$
$\displaystyle=\sqrt{\frac{2k}{e^{kL}-1}}\frac{k\alpha_{n}}{\sqrt{(E_{n}+q)^{2}+(k\alpha_{n})^{2}}}e^{\frac{k|y|}{2}}\sin\left[\alpha_{n}\left(1-e^{k|y|}\right)\right]\mathrm{sign}(y)\leavevmode\nobreak\
,$ (72) $\displaystyle\phi_{n}^{\downarrow}(y)$
$\displaystyle=\sqrt{\frac{2k}{e^{kL}-1}}\frac{E_{n}+q}{\sqrt{(E_{n}+q)^{2}+(k\alpha_{n})^{2}}}e^{\frac{k|y|}{2}}\cos\left[\alpha_{n}\left(1-e^{k|y|}\right)\right]\leavevmode\nobreak\
,$ (73)

are obtained for $\eta=-1$, and where the remaining constant was set by the
normalization of the wavefunction

$\int_{0}^{L}dy\tilde{\phi}_{n}(q,y)^{\dagger}\tilde{\phi}_{n}(q,y)=1\leavevmode\nobreak\
.$ (74)

The solution for the particular case of $n=0$ has only a lower component, and
is given by

$\phi_{0}^{\downarrow}(y)=\sqrt{\frac{k}{e^{kL}-1}}e^{\frac{k|y|}{2}}\mathrm{sign}(E_{n}+q)\leavevmode\nobreak\
.$ (75)

### QW explicit time step

Making use of the equations that define the QW, Eqs. (37, 38, 39) and (40),
one can recast the evolution of $|\chi_{j}\rangle$ as a recurrence relation
relating the spinor components Eq. (41) at two consecutive time steps. We
arrive at

$\displaystyle\chi_{j+1,r,s}^{\uparrow}=$
$\displaystyle-\frac{i}{2}e^{i\theta(y)}\left[s\left(y+\frac{\epsilon}{2}\right)+s\left(y-\frac{\epsilon}{2}\right)\right]\chi_{j,r+1,s}^{\uparrow}-\frac{1}{2}\left[s\left(y+\frac{\epsilon}{2}\right)-s\left(y-\frac{\epsilon}{2}\right)\right]\chi_{j,r-1,s}^{\downarrow}$
$\displaystyle+\frac{1}{2}f(y)f(y+\epsilon)c\left(y+\frac{\epsilon}{2}\right)\chi_{j,r+1,s+1}^{\uparrow}+\frac{1}{2}f(y)f(y-\epsilon)c\left(y-\frac{\epsilon}{2}\right)\chi_{j,r+1,s-1}^{\uparrow}$
$\displaystyle+\frac{i}{2}f(y)f^{*}(y+\epsilon)c\left(y+\frac{\epsilon}{2}\right)\chi_{j,r-1,s+1}^{\downarrow}-\frac{i}{2}f(y)f^{*}(y-\epsilon)c\left(y-\frac{\epsilon}{2}\right)\chi_{j,r-1,s-1}^{\downarrow}\leavevmode\nobreak\
,$ (76)

for the upper component, where we recall that $y=\epsilon s$, and we defined
$e^{\pm i\theta(y)}=c(y)\pm is(y)$. For the lower component one finds

$\displaystyle\chi_{j+1,r,s}^{\downarrow}=$
$\displaystyle\frac{i}{2}e^{-i\theta(y)}\left[s\left(y+\frac{\epsilon}{2}\right)+s\left(y-\frac{\epsilon}{2}\right)\right]\chi_{j,r-1,s}^{\downarrow}-\frac{1}{2}\left[s\left(y+\frac{\epsilon}{2}\right)-s\left(y-\frac{\epsilon}{2}\right)\right]\chi_{j,r+1,s}^{\uparrow}$
$\displaystyle+\frac{1}{2}f^{*}(y)f^{*}(y+\epsilon)c\left(y+\frac{\epsilon}{2}\right)\chi_{j,r-1,s+1}^{\downarrow}+\frac{1}{2}f^{*}(y)f^{*}(y-\epsilon)c\left(y-\frac{\epsilon}{2}\right)\chi_{j,r-1,s-1}^{\downarrow}$
$\displaystyle-\frac{i}{2}f^{*}(y)f(y+\epsilon)c\left(y+\frac{\epsilon}{2}\right)\chi_{j,r+1,s+1}^{\uparrow}+\frac{i}{2}f^{*}(y)f(y-\epsilon)c\left(y-\frac{\epsilon}{2}\right)\chi_{j,r+1,s-1}^{\uparrow}\leavevmode\nobreak\
.$ (77)

We notice that the upper components are displaced in one direction along the
$x$ dimension, while the lower components are displaced in the opposite
direction.

### Mode decomposition of the freely propagating distribution

The stationary states found above form an orthonormal basis, in the continuum
limit, that allow for a decomposition of any function along the $y$
coordinate, for a given value of $q$. They can also be used, after a proprer
discretization, in the lattice on which the QW is defined. Following this
idea, we introduced the decomposition in Eq. (46), which is a function in the
space of $q$, the lattice quasimomentum along the $x$ coordinate. For this
quasimomentum space, the spinor components are related to Eq. (41) via a
discrete Fourier transform

$\tilde{\chi}_{j,s}(q)=\sum_{r}e^{-iq\epsilon
r}\chi_{j,r,s}\leavevmode\nobreak\ .$ (78)

Making use of

$\sum_{r}e^{ix(q-q^{\prime})}=\frac{2\pi}{\epsilon}\delta(q-q^{\prime})\leavevmode\nobreak\
,$ (79)

and the orthonormality condition (74) on the grid

$\epsilon\sum_{s}\tilde{\phi}_{n}(q,\epsilon
s)^{\dagger}\tilde{\phi}_{m}(q,\epsilon s)=\delta_{n,m}\leavevmode\nobreak\ ,$
(80)

the coefficients can be obtained as

$\beta_{n}(q,t)=\epsilon^{2}\sum_{s}\tilde{\phi}_{n}(q,\epsilon
s)\tilde{\chi}_{j,s}(q)\leavevmode\nobreak\ .$ (81)

The coefficients of the freely propagating distribution with $x=t$ are

$\beta_{n}(q,t)=\epsilon^{2}\sum_{s}\tilde{\phi}_{n}(q,\epsilon
s)e^{-iqt}\chi_{j,j,s}\leavevmode\nobreak\ ,$ (82)

while, for $x=-t$, they read as

$\beta_{n}(q,t)=\epsilon^{2}\sum_{s}\tilde{\phi}_{n}(q,\epsilon
s)e^{iqt}\chi_{j,-j,s}\leavevmode\nobreak\ .$ (83)

From the normalization condition of the spinor on the grid

$\epsilon^{2}\sum_{r,s}\chi_{j,r,s}^{\dagger}\chi_{j,r,s}=1\leavevmode\nobreak\
,$ (84)

and making use of the definition (46), it can be shown that the mode
coefficients satisfy

$\sum_{n}\int_{-\pi/\epsilon}^{\pi/\epsilon}\frac{dq}{2\pi}\left|\beta_{n}(q,t)\right|^{2}=1\leavevmode\nobreak\
,$ (85)

which can be expressed in terms of the integrated coefficients (47) as

$\sum_{n}B_{n}(t)=1\leavevmode\nobreak\ .$ (86)

### High $kL$ limit of the QW time step and limiting entropy

In the limit of a high warp factor $kL$, the exponential $e^{-A(L)}$ becomes
very small, so that the QW discrete time recursive evolution Eqs. (76,77) can
be expanded up to the lowest order in this factor, giving

$\displaystyle\chi_{j+1,r,s}^{\uparrow}=$
$\displaystyle\chi_{j,r+1,s}^{\uparrow}\leavevmode\nobreak\ ,$
$\displaystyle\chi_{j+1,r,s}^{\downarrow}=$
$\displaystyle\chi_{j,r-1,s}^{\downarrow}\leavevmode\nobreak\ .$ (87)

Although this expansion is only valid for values of $y$ close to $L$, it is
still accurate enough for the initial condition located at $y=L/2$. As
discussed in the main text, the asymptotic value of the entanglement entropy
decreases as $kL$ is increased. Therefore, the minimum value of the entropy is
reached in the limit $e^{-A(L)}\approx 0$. The initial condition
$\chi_{0,r,s}=\delta_{r,0}\delta_{s,s_{0}}C_{0}$ can be iterated with the help
of Eqs. (87) to produce the explicit time evolution

$\displaystyle\chi_{j,r,s}^{\uparrow}=$ $\displaystyle
C_{0}^{\uparrow}\delta_{r,j}\delta_{s,s_{0}}\leavevmode\nobreak\ ,$
$\displaystyle\chi_{j,r,s}^{\downarrow}=$ $\displaystyle
C_{0}^{\downarrow}\delta_{r,-j}\delta_{s,s_{0}}\leavevmode\nobreak\ .$ (88)

The corresponding reduced density matrix becomes time-independent and
diagonal:

$\rho_{c}(t)=diag(|C_{0}^{\uparrow}|^{2},|C_{0}^{\downarrow}|^{2})$ (89)

from which the minimum value of the entropy can finally be obtained:

$S_{\mathrm{min}}=-|C_{0}^{\uparrow}|^{2}\log_{2}|C_{0}^{\uparrow}|^{2}-|C_{0}^{\downarrow}|^{2}\log_{2}|C_{0}^{\downarrow}|^{2}\leavevmode\nobreak\
.$ (90)