Emmanuel Baccelli

1 Pervasive Monitoring is a Myth

Pervasive monitoring does not exist for one good reason: if monitoring does happen, as revealed by recent whistleblowing that highlighted the scale of NSA surveillance activities [START_REF] Greenwald | Edward Snowden: the whistleblower behind the NSA surveillance revelations[END_REF], it happens in very few key central locations in the network. Monitoring, as we know it so far, is thus by no means pervasive in the geographical or topological sense. One might further speculate that such monitoring will probably never become pervasive because it would simply be too expensive to put in place.

This observation leads to a fundamental question: how could NSA surveillance reach such an incredibly large scope, to the point where it seems pervasive, while in fact, they only monitor the network from a few select locations? The answer to this question is simple: the NSA's monitoring could scale because we are addicted to centralization.

From the physical point of view, we made ourselves absolutely dependent on the deployed, fixed infrastructure. Whenever we communicate, the quasi totality of the traffic relies upon going through a central access point of some sort. This approach to networking certainly offers substantial advantages (a good performance/cost compromise for starters), but on the other hand it offers cheap and efficient vantage points for whoever seeks large scale monitoring opportunities. It does not have to be so.

Similarly, from the logical point of view, we make ourselves more and more dependent on cloud services. Whenever we store, or access our data, we rely on a centralized architecture that concentrates everything in data centers. Again, this approach offers great advantages (e.g. easier access from anywhere, lower cost reliability and efficiency), but on the other hand it offers more cheap and efficient vantage points for would-be spies. It does not have to be so either.

Making Massive Surveillance Harder

At this point, it becomes apparent that there are essentially three ways to make large scale systematic surveillance a much harder task.

The first type of solution is to significantly increase the physical security of potential vantage points. While this type of solution is obviously necessary, history shows that there are always cracks to infiltrate if the will to break in is there. And it will be there. So this type of solution is by no means sufficient.

The second type of solution is to significantly increase the default levels of encryption used over the network. For starters, the default behavior should be to encrypt "door-to-door" and don't let communication be decrypted by intermediate devices. Furthermore the default cryptographic mechanisms should be significantly harder to decipher. Note that to ensure this actually happens consistently, the use of open source software should be mandatory concerning such cryptographic mechanisms. On one hand, as highlighted in [START_REF] Schneier | Security Risks of Embedded Systems[END_REF], open-source bundled with automatic updates can slow down or stop obsolescence altogether, and on the other hand, open source is best known guarantee against malware and potential backdoors. However, if one of the physical end points of the communication is a vantage point of potentially high interest (say a data center), we quickly come back to the problems of the first type of solution, i.e. there are always cracks to infiltrate if the will to break in is there. And it will be there, with an extremely intense focus, because this is a potential vantage point -and with an intense focus one can do wonders.

The third type of solution is of an entirely different nature, which is the main purpose of this document. Instead of hardening the current architecture and increasing the security of high profile targets, another category of approach could be to aim for target dispersal, as suggested recently by B. Schneier in [START_REF] Schneier | What we know and what we do not know[END_REF]. Target dispersal would eliminate "default" vantage points and thus naturally disable systematic mass surveillance. De facto, surveillance efforts would be forced to be more specific and personalized, and thus more directly accountable for.

To that end, stateless networking approaches should be developed and employed, i.e. approaches that do not rely on central entities (infrastructure-based), but rather on spontaneous interaction between autonomous peer entities, as locally as possible (topologically and/or geographically). The network architecture should evolve towards a mode of operation where all possible stateless solutions are tried first, before considering any infrastructure-based approach -again, the goal being to increase target dispersal. In fact, several examples of stateless approaches have appeared over the last decade, at different layers. These include for instance peer-to-peer (P2P) networks or WebRTC [START_REF]Real-Time Communication in WEB-browsers IETF Working Group[END_REF] at the application layer, and multi-hop spontaneous wireless networks at the network layer [START_REF] Cordero | Enabling Multihop Communication in Spontaneous Wireless Networks[END_REF], such as mobile ad hoc networks [START_REF] Corson | Mobile Ad hoc Networking (MANET): Routing Protocol Performance Issues and Evaluation Considerations[END_REF], wireless sensor networks, vehicular networks, or wireless mesh networks.

While P2P networks have been massively deployed and adopted (WebRTC will likely enjoy the same fate) spontaneous wireless networks have not yet been widely adopted. There has however been quite some work in the research community over the past decade on this topic, which has identified issues that must (still) be tackled in order to efficiently adapt IP protocols for these new networks. Actually, the community has already produced a number of new protocol specifications dedicated to spontaneous wireless networking operation, such as [START_REF] Clausen | Optimized Link State Routing Protocol (OLSR)[END_REF] [6] [9], [START_REF] Winter | RPL: IPv6 Routing Protocol for Low-Power and Lossy Networks[END_REF] among others. Some of these protocols have even been deployed in specific contexts. Several european and american military applications have deployed using mobile ad hoc routing to power spontaneous, on-the-move local communications between elements on-site. In another context, several cities in Europe and in the USA have deployed wireless community mesh networks (e.g. [18] [19]), which use these protocols. Last, but not least, the highly anticipated Internet of Things (IoT) is heavily based on spontaneous wireless networking and the outcome of efforts such as 6LoWPAN [START_REF]IPv6 over Low power WPAN (6lowpan) IETF Working Group[END_REF] or ROLL [START_REF]Routing Over Low power and Lossy networks (roll) IETF Working Group[END_REF].

Indeed, merging spontaneous wireless networks with traditional, operated networks presents a number of advantages, aside of target dispersal to counter massive surveillance. First, it could offer operator infrastructure offloading [START_REF] Baccelli | Proceedings of the 3rd MANIAC Challenge[END_REF]. Second, the resulting network would systematically maximize connectivity at marginal cost. Third, it offers natural coverage extension for the network access infrastructure already deployed. And last, it offers increased resilience of the network in face of infrastructure outage.

There are of course lingering issues that do not allow off-the-shelf, spontaneous wireless networking at large scale. Among others, such issues include (i) the absence of efficient standard router auto-configuration schemes in this context [START_REF] Baccelli | IP Addressing Model in Ad Hoc Networks[END_REF], which has profound implications, (ii) the lack of dedicated, optimized link layer technology (which has focused on infrastructure based networking for decades), (iii) no standard key exchange protocol or alternatives for efficient authentication method to date, and (iv) some fundamental differences in terms of characteristics compared to traditional networks, such as throughput capacity vs number of nodes in the network.

However, in light of the recent events revealing massive and systematic surveillance [START_REF] Greenwald | Edward Snowden: the whistleblower behind the NSA surveillance revelations[END_REF], there is a strong argument to decentralize our networking paradigms, and in this realm, spontaneous wireless networking at the networks layer is an asset that should not be overlooked. 3 Position: More Spontaneity Please

The position highlighted in this paper is that in order to maximize target dispersal to counter systematic massive surveillance, we should not overlook the full potential of decentralized network paradigms. This full potential includes not only mantras such as "think hard before you use the cloud ", or "let's bypass the service provider " but also alternative techniques at the network layer, such as multi-hop spontaneous wireless networking, which could significantly reduce our dependence on the infrastructure -which will remain, whether we want it or not, a desirable vantage point for mass monitoring. There are a number of issues that need to be resolved or alleviated towards native and massive integration of spontaneous wireless networks in the currently deployed IP architecture and infrastructure-based networks (massive IPv6 adoption could already help somewhat). However, the gains obtained in terms of target dispersal alone could be worth this cost -and there are other gains too.

  Contrary to infrastructure-based approaches, which are prone to monitoring, spontaneous wireless networking uses communication links that are local and volatile, i.e. unless one is physically present at the time and location of the communication, one must abandon all hope of monitoring anything. While the wireless nature of the communication may on the other hand facilitate eavesdropping, the fact that one has to be there at the time/place of the communication significantly hampers mass surveillance, and cryptographic techniques can provide a privacy equivalent to what is achieved on wire. The honest, technical question that we should ask ourselves at this point is: are the lingering issues of spontaneous wireless networking integration in the IP stack worth solving/improving or not, taking into account our goal to maximize surveillance target dispersal?


Emmanuel Baccelli

Oliver Hahm

Mesut Günes

Matthias Wählisch

Thomas C Schmidt

RIOT OS: Towards an OS for the Internet of Things

The Internet of Things (IoT) is characterized by heterogeneous devices. They range from very lightweight sensors powered by 8-bit microcontrollers (MCUs) to devices equipped with more powerful, but energy-efficient 32-bit processors. Neither traditional operating systems (OS) currently running on Internet hosts, nor a typical OS for sensor networks are capable to fulfill the diverse requirements of such a wide range of devices. To leverage the IoT, redundant development should be avoided and maintenance costs should be reduced. In this paper we revisit the requirements for an OS in the IoT. We introduce RIOT OS, an OS that explicitly considers devices with minimal resources but eases development across a wide range of devices. RIOT OS allows for standard C and C++ programming, provides multithreading as well as real-time capabilities, and needs only a minimum of 1.5 kB of RAM.

I. INTRODUCTION

Billions of heterogeneous devices such as sensor nodes, home appliances, smartphones, and vehicles are expected to be interconnected by spontaneous wireless networks or power line communication, thus giving birth to the IoT [START_REF] Ashton | That 'Internet of Things' Thing[END_REF]. Such devices, though extremely constrained in terms of computing power, available memory, communication, and energy capacities, are however expected to fulfill the requirements of cyperphysical systems: (i) reliability, (ii) real-time behavior, and (iii) an adaptive communication stack to integrate the Internet seamlessly.

An OS for IoT devices should fulfill these requirements and run on a wide spectrum of hardware, ranging from nodes based on low-power MCUs, to nodes powered by new generations of energy-efficient 32-bit processors. None of the existing OS -neither a lightweight OS targeting Wireles Sensor Networks (WSNs) nor a full-fledged OS -is capable to fulfill requirements so diverse. Hence, in order to avoid redundant developments and maintenance costs, a new, unifying type of OS is needed, which is the subject of this poster.

Problem Statement: Ideally, the capabilities of a fullfledged OS should be available on all IoT devices. The following will focus on the characteristics of Linux, since it is a good example for an open source OS among the major fullfledged OS. Compared to a typical lightweight OS targeting WSNs, Linux is more developer-friendly: numerous available system libraries, network protocols or algorithms, and nearzero learning curve in the sense that developers can code in standard C or C++. However, Linux's minimal requirements in terms of CPU and memory do not fit constrained IoT devices powered by small MCUs. While efforts have attempted at pushing down these requirements [START_REF] Mccullough | uCLinux for Linux Programmers[END_REF], we argue that Linux has not been designed for the IoT and cannot fulfill strict energy efficiency. For these reasons Linux cannot be expected to become the one OS to rule them all in the IoT.

On the other hand, the trade-offs that enable a typical lightweight OS targeting WSNs to run on the most constrained IoT devices make it significantly less developer-friendly and its use on more powerful IoT devices will result in a less energyefficient implementation, while not exploiting devices' full capabilities. The dominant WSN OS, Contiki and TinyOS [START_REF] Farooq | Operating systems for wireless sensor networks: A survey[END_REF], follow an event driven design, which is useful for typical WSN scenarios, but exhibit drawbacks for efficient and functional networking implementations. For example, in a typical WSN scenario it is sufficient to process the tasks sequentially in a loop, but doing so limits the windows sizes to at most one TCP connection, due to constrained memory. An alternative system architecture is desirable, providing capabilities of a modern, full-fledged OS, such as native multi-threading, hardware abstraction, dynamic memory management. However, these types of systems have so far been considered too complex for IoT devices.

II. RIOT: AN OS FOR SMALLER THINGS IN THE INTERNET

Design Aspects for an IoT OS: Fundamentally, an OS is characterized by the following key design aspects: the structure of the kernel, the scheduler, and the programming model. The kernel can either (i) be built in a monolithic fashion, (ii) follow a layered approach, or (iii) implement the microkernel architecture. The choice of the scheduling strategy is tightly bound to real-time support (or the lack thereof), the support for different task priorities, or the supported degree of user interaction. Finally, the programming model defines whether (i) all tasks are executed within the same context and have no segmentation of the memory address space, or (ii) every process can run in its own thread and has its own memory stack. The programming model is also linked to the available programming languages for application developers.

Comparing Existing Solutions: Based on these observations, we compare Contiki, Tiny OS, and Linux, and we derive first principles for an OS for the IoT. Contiki follows a modular concept close to the layered approach, while Tiny OS consists of a monolithic kernel, as Linux. The scheduling in Contiki is purely event driven, similar to that in TinyOS, where a FIFO strategy is used. Linux on the other hand, uses a scheduler, which guarantees a fair distribution of processing time. The programming models in Contiki and TinyOS are based on the event driven model, in a way that all tasks are executed within the same context, although they offer partial multithreading support. Contiki uses a subset of the C programing language, where some keywords cannot be used, while TinyOS is written in a C dialect called nesC. Linux, on the other hand, supports real multi-threading, is written in standard C, and offers support for a wide range of different programming and scripting languages. Because of these design trade-offs, TinyOS and Contiki are thus lacking several key developerfriendly features: standard C and C++ programability, standard multi-threading, and real-time support (see Table I). An OS leveraging a successful, large-scale deployment of IoT devices should support these functions. From the developer perspective, this means a powerful, hardware-independent API. RIOT Architecture Overview: RIOT OS aims at bridging the gap we observed between OS for WSNs and traditional full-fledged OS currently running on Internet hosts. It is based on design objectives including energy-efficiency, small memory footprint, modularity, and uniform API access, independent of the underlying hardware.

OS Min RAM Min ROM C Support C++ Support Multi-Threading MCU w/o MMU Modularity Real-Time Contiki <2kB <30kB • ✘ • ✓ • • Tiny OS <1kB <4kB ✘ ✘ • ✓ ✘ ✘ Linux ~1MB ~1MB ✓ ✓ ✓ ✘ • • RIOT ~1.5kB ~5kB ✓ ✓ ✓ ✓ ✓ ✓
RIOT implements a microkernel architecture inherited from FireKernel [START_REF] Will | A real-time kernel for wireless sensor networks employed in rescue scenarios[END_REF], thus supporting multi-threading with standard API. In addition to FireKernel's original features, RIOT adds support for C++ -enabling powerful libraries such as the Wiselib algorithm framework -and provides a TCP/IP network stack. Advantages of the RIOT architecture thus include: (i) high reliability and (ii) a developer-friendly API. The modular microkernel structure of RIOT makes it robust against bugs in single components. Failures in the device driver or the file system, for example, will not harm the whole system. RIOT allows developers to create as many threads as needed and distributed systems can be easily implemented by using the kernel message API. The amount of threads is only limited by the available memory and stack size for each thread, while the computational and memory overhead is minimal.

RIOT Technical Details: To fulfill strong real-time requirements RIOT enforces constant periods for kernel tasks (e.g., scheduler run, inter-process communication, timer operations). An important prerequisite for guaranteed runtimes of O(1) is the exclusive use of static memory allocation in the kernel. Yet, dynamic memory management is provided for applications. We achieve constant runtime of the scheduler by using a fixed-sized circular linked list of threads. Constant runtime of the timer operations is obtained by exploiting the fact that MCUs typically provide multiple compare registers.

It is mandatory to maximize the duration spent in deep-sleep modes, in order to implement energy-efficiency also for more powerful IoT devices, RIOT thus introduces a scheduler that works without any periodic events. Whenever there are no pending tasks, RIOT will switch to the idle thread, which determines the deepest possible sleep mode depending on peripheral devices in use. Only interrupts (external or kernelgenerated) wake up the system from idle state. Low complexity of kernel functions is a main factor for the energy efficiency of an OS. Therefore, the duration and occurrence of context switching have to be minimized. In RIOT context switching is performed in two cases: (i) a corresponding kernel operation itself is called, e.g., a mutex locking or creation of a new thread, or (ii) an interrupt causes a thread switch. The first case will occur rarely. For example, every thread is usually created once. Hence, it is important to reduce the processing time in case of a thread switch. Therefore, RIOT's kernel provides a minimized scheduler, when it gets called out of an interrupt service routine. In that case, saving the old thread's context is not required and thus a task switch can be performed in very few clock cycles.

III. AVAILABLE CODE & FUTURE WORK

Despite this sophisticated architecture and efficient scheduling, RIOT has a low memory footprint. Available open source RIOT code [START_REF]RIOT Operating System[END_REF] requires less than 5 kByte of ROM and less than 1.5 kByte of RAM for a basic application on MSP430, for instance. RIOT uses a multi-threaded programming model in combination with standard ANSI C code and a common POSIX-like API for all supported hardwarefrom 16-bit microcontrollers to 32-bit processors. Hence, for projects involving heterogeneous IoT hardware, it is possible to build the whole software system upon RIOT and easily adopt existing libraries. Moreover, the availability of several networking protocols including the latest standards of the IETF for connecting constrained systems to the Internet (6LoWPAN, RPL) make RIOT IoT-ready. On-going work includes full POSIX compliance and porting to various IoT platforms.

TABLE I KEY

 I CHARACTERISTICS OF CONTIKI, TINYOS, LINUX, AND RIOT. (✓) FULL SUPPORT, (•) PARTIAL SUPPORT, (✘) NO SUPPORT. THETABLE COMPARES THE OS IN MINIMUM REQUIREMENTS IN TERMS OF RAM AND ROM USAGE FOR A BASIC APPLICATION, SUPPORT FOR PROGRAMMING LANGUAGES, MULTI-THREADING, MCUS WITHOUT MEMORY MANAGEMENT UNIT (MMU), MODULARITY, AND REAL-TIME BEHAVIOR.


Grasha Jacob
email: grasharanjit@gmail.com

Dr A Murugan

On the Construction of Doubly Even Order Magic Squares

Keywords: Square matrix, Magic Square, Doubly even order, Magic Constant

Magic squares have been an enthralling topic in mathematics for centuries. They are formed by filling in all the cells of a square matrix with the numbers starting from one so that the sum of all rows, columns, and diagonals is the same. Magic squares have applications in entertainment, music and even cryptography. This paper focuses on a simple and easy method to construct doubly even magic squares.

INTRODUCTION

A magic square of order n is a square matrix or array of n 2 numbers such that the sum of the elements of each row and column, as well as the main diagonals, is the same number, called the magic constant denoted by (M). Generally, the entries are thought of as the natural numbers 1, 2, ..., n 2 , where each number is used exactly once. Magic Squares exist for all orders, n ≥ 1 with an exception when n = 2. The case n = 1 is trivial as it consists of a single cell containing the number 1 alone. The constant sum in every row, column and diagonal is called the magic constant or magic sum. The magic constant of a normal magic square depends only on n. For normal magic squares of order n = 3, 4, 5, 6, 7, 8 … the magic constants are 15, 34, 65, 111, 175, 265… respectively. Magic squares are classified into three types: odd, doubly even and singly even. A doubly even magic square is a square matrix of order n, where n is divisible by four. A singly even magic square is a square matrix of order n, where n is even but not divisible by four.

PROPERTIES OF MAGIC SQUARES

 The magic constant of a normal magic square depends only on n, where n is the order of the magic square. The magic constant is given by the formula,

(M) = n(n 2 +1)/2.
For normal magic squares of order n = 3, 4, 5, 6, 7, 8 … the magic constants are 15, 34, 65, 111, 175, 265… respectively.  A magic square will remain magic if any number is added to every number of a magic square.  A magic square will remain magic if any number multiplies every number of a magic square.  The transpose of a magic square is also a magic square [START_REF] Stephens | Matrix Properties Of Magic Squares[END_REF].  A magic square will remain magic if two rows, or columns, equidistant from the center are interchanged.

 An even order magic square will remain magic if the diagonally opposite quadrants are interchanged.  An even order magic square will remain magic if the columns or the rows are exchanged. If the 1st and the 2nd columns/rows are exchanged and the 3rd and the 4th columns/rows are exchanged, this is called a 1-2-3-4 exchange.

A most-perfect (MP) magic square MM is of doubly-even order (n = 4, 8, . . .) and has the following properties:

 two elements that are n/2 elements apart along all diagonals (including broken ones in both directions) sum to n 2 +1.  the elements of all 2 by 2 subsquares (including broken top-bottom, broken leftright, and the four corners) sum to the constant 2n 2 +1.

MAGIC SQUARES OF DOUBLY EVEN ORDER

A doubly even magic square is constructed by drawing x's through each 4 x 4 subsquare, filling all squares in sequence and then replacing each entry aij on a crossedoff diagonal by n x n + 1 -a ij or, equivalently, reversing the order of the crossed-out entries. For n = 8, the crossed-out numbers are originally 1, 4, ... , 61, 64, so entry 1 is replaced with 64, 4 with 61, etc as in figure 1.

Figure 1. Magic Square of Order 8

PROPOSED METHOD

There are many ways to construct magic squares, but the standard way is to follow certain formulas for generating the patterns. In this paper, a doubly even magic square is constructed by drawing x's through the main diagonals of each 4 x 4 subsquare, filling all squares in sequence and then replacing each entry a ij other than the crossed-off diagonal by n x n + 1 -a ij and having crossed off entries as such. Figure 2 symbolizes the proposed method.

Figure 2. Construction of Magic Squares

Algorithm Magic Square

Input : Order of matrix, n Output : Magic Square, A 1. Read the value of n, where n is doubly even.

Declare a matrix A of order n X n. Fill all the cells of the matrix with zeros. 3. Divide the matrix into blocks of size 4 x 4. 4. Fill the diagonal of each 4 X 4 block with the corresponding cell value, c ij . 5. Trace through each row and fill all the cells of the entire matrix with zero entries

with the value of n x nc ij . End Magic Square

4.1

Steps in the construction of Doubly Even order Magic Square Step 1 Read the value of n, where n is doubly even. Let us assume n = 8.

Step 2 Consider a matrix A of order 8 X 8. Fill all the cells of the matrix with zeros. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Step 3 Divide the matrix into blocks of size 4 x 4.

Step 4 Fill the diagonal of each 4 X 4 block with the corresponding cell value, c ij .

Note:

i) For the main diagonal of a 4 x 4 block, i=j, where i and j represent the row and column of the matrix. For the diagonal in the opposite direction, i+j = 5.

ii) The cell value of each cell is obtained by the formula c ij =(i-1) * n + j, where, i,j represent the row and column of the matrix and n is the order of the matrix.

Step 5 Trace through each row and fill all the cells of the entire matrix with zero entries with the value of n x nc ij . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

4.2

Matlab Code for Magic Square of order n, where n is doubly even

CONCLUSION

In this paper, a simple and easy method of constructing a doubly even magic square is proposed and implemented. The algorithm and the code in Matlab are also given so that the applications of magic squares can be explored in detail.


Soo-Young Lee

Jérémy Le Deunff

Muhan Choi
email: mhchoi@ee.knu.ac.kr

Roland Ketzmerick

Quantum Goos-Hänchen shift and tunneling transmission at a curved step potential

Keywords: numbers: 03.65.Vf, 03.65.Xp, 42.25.Gy

We study the quantum Goos-Hänchen (GH) shift and the tunneling transmission at a curved step potential by investigating the time evolution of a wave packet. An initial wave packet is expanded in terms of the eigenmodes of a circular step potential. Its time evolution is then given by the interference of their simple eigenmode oscillations. We show that the GH shift along the step boundary can be explained by the energy-dependent phase loss upon reflection, which is defined by modifying the one-dimensional (1D) effective potential derived from the 2D circular system. We also demonstrate that the tunneling transmission of the wave packet is characterized by a freespace image distant from the boundary. The tunneling transmission exhibits a rather wide angle divergence and the direction of maximum tunneling is slightly rotated from the tangent at the incident point, which is consistent with the time delay of the tunneling wave packet computed in the 1D modified effective potential.

I. INTRODUCTION

When an optical beam is totally reflected at a dielectric interface, the reflected beam comes out at a point shifted along the interface from the incident point. This deviation is called Goos-Hänchen (GH) shift [1]. This arises due to the interference of its plane wave components whose phase loss upon the reflection varies with their incident angle. Although it has been first observed in optical beams, the GH shift (or effect) is a generic phenomenon of wave mechanics and can be found in various physical systems such as graphene [2,3], normal/superconductor interface [4], optical beam [5,6], dielectric microcavity [7][8][9] or in the dynamics of neutron [10].

Recently, the GH effect has been studied with electron waves in condensed matters, p/n-doped interface in graphene [2] and normal/superconductor (NS) interface [4]. Note that these works have investigated the GH effect appearing at a planar interface, where it is not difficult to find the phase loss, occurring upon reflection, for a plane wave component. However, if we deal with a curved interface, the plane wave decomposition becomes useless as the plane wave components are no longer satisfying the matching conditions at the curved boundary. Thus a straightforward extension of the method used in the planar interface cannot be applied to understand the GH shift appearing along the curved boundary and it becomes necessary to provide a qualitative/quantitative description in such a case. The GH shift at a curved interface has been studied by a few authors and only for optical systems. The contribution of the GH shift in a resonance mode of a circular microcavity has been esti-mated by comparing the decay rate expected from a multiple reflection with the eigenvalue of a resonance mode [8]. Also, the deviation of an optical-beam path from the specular ray reflection has been discussed when an incident beam comes from outside the cavity and reaches the circular interface [9].

In addition to the GH effect, the curvature of the boundary induces tunneling transmission which cannot be observed with the simple planar interface. It will occur at a range of energy just higher than the height of the potential and corresponds to dynamical tunneling [11], even though it can be interpreted as tunneling through the barrier of an effective 1D potential, as we shall see later on. It is noted that both the GH shift and the tunneling arise from the intrinsic nature of (quantum) waves and have no classical correspondence. Recently, considerable efforts have been made to understand tunneling phenomena based on the associated classical dynamics in various systems such as quantum map [12][13][14][15], optical microcavities [16][17][18][19] or quantum potentials [20,21].

In this paper, we study the time evolution of wave packets in the 2D circular step potential. We investigate the quantum GH shift and the tunneling transmission observed in the wave packet dynamics. Using the eigenmodes (bound modes and resonance modes) of the system introduced in Section II, we construct an initial wave packet located deep inside the circular cavity, and its subsequent time evolution is then described by the interference of simple oscillations in time of individual eigenmodes (Sec. III). Section IV is dedicated to the GH shift. We evaluate the phase change upon reflection by modifying partly the effective potential derived from the 2D system. It follows from the assumption that the GH shift is relevant only to the barrier region of the potential. It is shown that the variation of the phase loss with respect to the energy is proportional to a time delay along the radial direction and is found to be consistent with the GH shift observed numerically. In Section V, we demonstrate that the transmitted (tunneling) wave packet exhibits a free-space image with a distance ∆ from the circular boundary. The distance ∆ can be understood as an average over individual distances ∆ j given by the angular momentum conservation in the j-th resonance mode. Then, we shall provide a qualitative explanation of the transmitted angle α T of the maximum tunneling direction in terms of the time delay of a 1D wave packet in a simplified 1D effective potential. Finally, we shall discuss the case of high-energy wave packets in Section VI.

II. EIGENMODES: BOUND AND RESONANCE MODES

We consider a 2D circular step potential with zero potential inside the circle and V 0 outside (see Fig. 2

(a)), i.e. V (r) = 0 if r < R ; V 0 if r ≥ R ,
where R is the radius of the circle. The rotational symmetry allows us to separate the time-independent Schrödinger equation

- 2 2m * ∇ 2 + V (r) Ψ(r, θ) = EΨ(r, θ) , (1)
into an angular and a radial equation by writing the wavefunctions Ψ(r, θ) as a product of the angular and radial solutions, respectively Θ(θ) and Φ(r). The solutions Θ(θ) of the angular part can be expressed as

Θ(θ) ∼ cos mθ, sin mθ, ( 2
)
where m is a natural integer corresponding to the angular momentum quantum number. The radial equation is written as

- 2 2m * ∂ 2 ∂r 2 + 1 r ∂ ∂r Φ(r) + V eff (r)Φ(r) = EΦ(r), (3)
where m * is the particle mass, and the effective 1D potential is defined as

V eff (r) = 2 m 2 2m * r 2 if r < R ; V 0 + 2 m 2 2m * r 2 if r ≥ R . (4)
For the sake of clarity, we will take m * = = 1, and R = 2 and V 0 = 5000 for the numerical calculations throughout the paper, unless other values are explicitly specified. An example of the effective potential V eff (r) with m = 120 is shown in Fig. 1 (a). It is easy to see that, in addition to the low energy range (E < V 0 ) associated with the bound modes, there is an energy range in which tunneling transmission can take place, i.e.

V 0 < E < V 0 + 2 m 2 /2m * R 2 .
The radial equation is nothing but the Bessel equation and the solution inside the circle is the Bessel function

Φ I (r) ∼ J m (kr) if r ≤ R , (5)
where k ≡ (2m * / 2 )E. Outside the circular cavity, the solution must be a decaying function when the energy E is less than the step height V 0 , while it would be a propagating function for E > V 0 . These solutions correspond to the modified Bessel function and the Hankel function, respectively

Φ II (r) ∼ K m (κr) if r ≥ R, E < V 0 ; H (1) m (k out r) if r ≥ R, E > V 0 , (6)
where κ ≡

(2m * / 2 )(V 0 -E) and k out ≡ (2m * / 2 )(E -V 0 ).
When E > V 0 we impose the outgoing-wave boundary condition by taking only the Hankel function of the first kind H m (k out r). Therefore, the resonance modes shall decay in time and thus are characterized by the decay rate γ corresponding to the imaginary part of the complex energy, E = ω -iγ/2. The complex energy indicates that we are dealing with a non-Hermitian system where the eigenmodes are not orthogonal to each other, as we shall discuss in the next section.

The solutions k mn or E mn (= 2 k 2 mn /2m * ) are extracted by imposing the boundary conditions,

Φ I (R) = Φ II (R), dΦ I dr (R) = dΦ II dr (R) .
Since our system is integrable, we have two good quantum numbers (m, n) specifying all the eigenmodes, where m and n are the angular momentum and radial quantum numbers, respectively. The k mn 's are shown in the (m, Re(k)) plane in Fig. 1 (b). The blue line denotes k B = m/R corresponding to the bottom of the effective potential. The red line is k V = (2m * / 2 )V 0 and separates the bound and resonance modes. Among the resonance modes, the modes below the green line,

k T = (2m * / 2 )V 0 + m 2 /R

III. INNER PRODUCT AND EIGENMODE EXPANSION

It is well known that the bound modes satisfy the orthogonality relation,

φi | φi ′ = δ ii ′ , (7)
where the tilde indicates that the mode is normalized as follows

| φi = 1 φ i |φ i |φ i . ( 8
)
The inner product is defined in the position representation as

φ i |φ i = d x φ * i ( x)φ i ( x) , (9)
where the star stands for the complex conjugate of the wavefunctions and x is a vector in the coordinate space, the (x, y)-plane in our case. However, the resonance modes |φ j are neither square integrable nor orthogonal, but hold the biorthogonality relation between the left ψj | and the right eigenmodes

| φj , ψj | φj ′ = δ jj ′ , (10)
where

| φj = 1 ψ j |φ j |φ j ; ( 11
) ψj | = 1 ψ j |φ j ψ j | . ( 12
)
We shall consider systems with time-reversal symmetry where the left eigenmode ψ j | is merely the complex conjugate of the right eigenmode [22], and the inner product is thus written as

ψ j |φ j = φ * j |φ j = d x φ 2 j ( x) . (13)
For the circular case, we can derive analytic expressions for the above inner products. With

Z m (x) ≡ J m (x) or H (1) m (x), one can use the indefinite-integral relation dx x Z 2 m (x) = 1 2 x 2 Z 2 m (x) -Z m-1 (x)Z m+1 (x) ,
to obtain the inner product of the resonance modes

φ * j |φ j = πR 2 2 J 2 m (z) -J m-1 (z)J m+1 (z) (14) - c 2 πR 2 2 H (1)2 m (z ′ ) -H (1) m-1 (z ′ )H (1) m+1 (z ′ ) ,
where the complex arguments are z ≡ kR and z ′ ≡ k out R, and c ≡ H

m (z ′ )/J m (z). A similar expression can be derived for the bound modes, replacing H

m (x) by K m (x),

φ i |φ i = πR 2 2 J 2 m (z) -J m-1 (z)J m+1 (z) (15) - c ′2 πR 2 2 K 2 m (z ′ ) -K m-1 (z ′ )K m+1 (z ′ ) ,
where z ≡ kR and z ′ ≡ κR, and c ≡ K m (z ′ )/J m (z). Note that k, κ, and c are real in the latter case.

Assuming the completeness relation

i | φi φi | + j | ψj φj | = 1 , (16)
where i and j run over the bound and resonance modes, respectively, an arbitrary state |Ψ can be expanded in terms of the eigenmodes of our system as follows

|Ψ = i φi |Ψ | φi + j ψj |Ψ | φj . ( 17
)
We consider now a 2D Gaussian wave packet with an average momentum k 0 ,

Ψ G ( x, t 0 ) = σ/π 1 + iσt 0 e -(σ/2)| x-x C | 2 +i k 0 •( x-x C )-iE 0 t 0 / 1+iσt 0 , ( 18
) where E 0 = | k 0 | 2 /2m * .
This is an exact solution of the free-particle time-dependent Schrödinger equation. The numerical calculations shall be carried out for σ = 100 and the position of the wave packet x C = (0, R) at t 0 = 0. We choose as an initial wave packet the Gaussian wave packet at t 0 = -1/k 0 , located deep inside the cavity, a unit length apart from the x C = (0, R) (see Fig. 2 (a)). In order to avoid any ambiguities in the determination of their positions and velocities, the parameters of all the considered wave packets will always be chosen such that their typical sizes are much smaller than the radius R of the potential. This initial wave packet would hit the point x C = (0, R) at the boundary after a time interval s 0 = 1/k 0 . For the following, we will conveniently set the time unit as s 0 .

Assuming a Gaussian wave packet with an angular momentum m 0 and a central wave number k 0 , as indicated by the large solid dots in Fig. 1 (b), then its central energy is 2 k 2 0 /2m * and the incident angle χ (see Fig. 2 (a)) is determined by the semiclassical relation sin χ = m 0 /k 0 R. The mean position of the initial wave packet Ψ

I ( x) = Ψ G ( x, t 0 = -1/k 0 ) is given by (x 0 , y 0 ) = (-m 0 /k 0 R, R -1 -(m 0 /k 0 R) 2 ).
The initial state can be expanded as

|Ψ I = i a i | φi + j b j | φj , (19)
where the overlap of the initial wave packet with the bound and resonance modes are respectively given by ), and presents an excellent agreement with the gaussian wave packet (Eq. ( 18)) as shown in Fig. 2 (b). It turns out that the reconstruction of the initial wave packet by Eq. ( 19) works well for the low-energy wave packets (A, B), but presents some discrepancies for the high-energy wave packets (C,D), as shown in Fig. 9 (a) and (b) where some modulation of the incident wave packets is seen.

IV. GH SHIFT IN WAVE PACKET DYNAMICS A. Numerical results

The time evolution of the initial wave packet is obtained by integrating the time-dependent Schrödinger equation and leads straightforwardly to

|Ψ(t) = i a i e -iEit/ | φi + j b j e -iEjt/ | φj . (20)
The wave packet would hit the boundary point (0, R) at t = 1 in the time unit of s 0 = 1/k 0 . The previous expression is nothing but Eq. ( 19) with the dynamical phase factors containing the eigenvalues E i and E j .

As an example, let us first consider the wave packet A with (m 0 , k 0 ) = (75, 75), which has an incident angle χ = π/6 and an energy E 0 = ( k 0 ) 2 /2m * = 2812.5 far below V 0 = 5000. We expand the initial wave packet in terms of 2474 eigenmodes (2428 bound modes and 66 resonance modes), which corresponds to the modes with an overlap |a i | (or |b i |) greater than 0.0005. The wave packet A at t = 0.5 and t = 1.5 is shown in Fig. 3 (a). The green dotted arrows indicate the path of the top of the wave packet and the GH shift is already clearly discernible. It is numerically evaluated at l GH ≃ 0.015.

For the wave packet B with (m 0 , k 0 ) = (120, 90) and χ ≃ 0.73 or 41.8 degrees, the same picture is shown in Fig. 3 (b). In this case, the GH shift is l GH ≃ 0.024. Although the energy E 0 = 4050 is less than V 0 , it contains many tunneling modes (608 tunneling modes) in addition to the bound modes (1374 bound modes) as shown in Fig. 1 (b), yielding an appreciable tunneling transmission that we shall discuss later on.

The shift l GH is computed using a more systematic evaluation of the wave packet position, by calculating the average position for a given t,

x av (t) = d x x |Ψ( x, t)| 2 d x |Ψ( x, t)| 2 , (21)
when the wave packet is sufficiently far away from the reflection point x C = (0, R). The wave packet positions at t = 0.4, 0.5, 0.6 and t = 1.4, 1.5, 1.6 are shown (blue solid dots) in Fig. 3 (c) and (d) for the wave packets A and B, respectively.

B. Theory for the GH shift at a curved boundary

In this section, we provide a simple explanation of the GH shift observed in the wave packet dynamics. In order to gain some insight into the theory, it is illuminating to mention first the simple case of the planar step potential. For a 1D step potential with a height V 0 , the reflection coefficient is expressed as F = e iφ when the energy E of the incident plane wave is less than V 0 . The phase of the reflection coefficient has an analytic expression

φ = 2 arctan E V 0 -E -π. ( 22
)
Then the wave packet dynamics exhibits a time delay upon reflection [4,24] which is given by the change of the phase φ with respect to E,

δt R = dφ dE = V 0 V 0 -E E + E V 0 -E , ( 23
)
and it is also known as the Wigner time scale [23]. Therefore, it is straightforward to see that the GH shift in the 22)) with V0 = 5000. The phases are nearly the same except for higher energy as shown in the inset.

2D planar case is characterized by a spatial shift l GH which is merely related to the time delay δt R by the expression

l GH = v y δt R , (24)
where v y is the velocity component parallel to the step boundary.

Returning to the circular case, let's take advantage of the 1D effective potential V eff (r) from which we want to extract the phase loss ϕ R under reflection in the same spirit as in the planar case. To do this, we assume that the time delay δt R is only determined by the shape of the barrier itself and the form of the potential where the incident wave packet evolves before hitting this barrier is not relevant. We can then modify the left part of potential into a constant potential so that we can consider plane wave solutions and the reflection coefficient F = |F |e iϕR in this region as illustrated in Fig. 4 (a). The reflection coefficient F is determined by matching the solutions at the frontier between the free region, i.e. the plane waves, and the barrier region, i.e. K m (κr) for E < V 0 and H (1) m (k out r) for E > V 0 . The figure 4 (b) shows the resulting phase loss ϕ R (solid red line) which is very similar to the phase loss φ of the step potential given by Eq. ( 22), except for a slight deviation in the high energy range of V 0 < E < V 0 + ( m) 2 /2m * R 2 (see the dashed green line in the inset). Then, similarly to the 2D planar case, defining the time delay in the radial direction as δt R = dϕR dE , the GH shift is given by

l GH = v θ δt R = ( m 0 R ) dϕ R dE , (25)
where v θ is the velocity of the wave packet projected along the θ-direction.

For the wave packet A with (m 0 , k 0 ) = (75, 75), Eq. ( 25) gives l GH = 0.0152 and t R ≃ 4.05 × 10 -4 ≃ 0.0304s 0 . The resulting trajectory is depicted by the red rectangles in Fig. 3 (c) where we used, for a better fit, a slightly increased reflection angle χ R = 1.01χ compared to the incident angle. This leads to an excellent agreement with the motion x av (t) of the wave packet computed numerically, illustrated by the blue dots on Fig. 3 (c). Similarly, for the wave packet B with (m 0 , k 0 ) = (120, 90) we get l GH = 0.0242 and t R ≃ 4.05 × 10 -4 ≃ 0.0363s 0 from Eq. ( 25). The reflection angle is also slightly modified χ R = 1.025χ and we see again a good agreement between the theory and the numerical results from the wave packet evolution as shown on Fig. 3 (d). The correction in the fit of the reflected angles can be justified as follows: on contrary to the 2D step potential case, in the circular cavity, the reconstructed wave packets (for instance A and B) contain tunneling-mode components which are responsible for a slight increase of the reflected angle. The reduction of smaller-incident-angle components, due to tunneling transmission, in the reflected wave packet results in the slight increase of the reflection angle χ R . In optics, this is called "Fresnel filtering effect" [25]. For a generic continuous 2D potential barrier, the GH shift would appear in a rather complex way due to the mixing with the change of dynamical phase of the rounded reflection at the continuous barrier.

V. TUNNELING TRANSMISSION A. Numerical results

As mentioned previously, the wave packet B contains many tunneling-mode components. So, we can expect to detect a transmitted wave packet outside of the circular region. In optical microcavities, the tunneling transmission has been characterized by its free-space image with a distance ∆ from the cavity surface [16,19,26]. The characterization by the distance ∆ is also pertinent in our case since both transmissions occur through the same mechanism, the tunneling process.

In order to determine the most probable direction of the tunneling transmission, we use the emission Husimi function used in Ref. [19], which is modified to be a wave-number independent function in the present case. We consider a projection line (h-axis) which has been rotated by an angle α and is a distance D apart from the origin as shown in Fig. 5

(a). The modified emission

Husimi function on the projection line is defined as

H E (D, h; α) = dh ′ Ψ(x(h ′ ), y(h ′ ))ξ(h, h ′ ) 2 ( 26
)
where Ψ(x, y) is the 2D wave function at a given time defined in Eq. ( 20) in position representation, h is the coordinate on the line and ξ(h, h ′ ) = It is worth noticing that the extrapolated position of the tunneling wave component at t = 1 is located at x ≃ 0 which is also the position of incident wave packet at t = 1 as illustrated in Fig. 5 (e). We can then estimate the time t * when the tunneling wave component is located at (x T , y T ) and we get t * ≃ 1 + 0.105 × (6.5/3.0) ≃ 1.227. Therefore, the tunneling wave component comes out with a time delay δt * ≃ 0.227s 0 at the position (x T , y T ) and a slightly rotated angle α T = -0.045.

B. Explanation of the distance ∆

In this section we aim to provide a semiclassical interpretation of the width ∆. We shall illustrate our purpose with help of the wave packet B. In the circular case, each individual resonance with (m j , k j ) exhibits its own distance ∆ j and has a clear physical meaning as it corresponds to the width of the effective potential Eq. ( 4) at the given energy E j = ( k j ) 2 /2m * . This is illustrated in Fig. 6 for the resonance mode (m j , k j ) = (120, 113.0 -i1.57 × 10 -6 ). This relation is obvious from the semiclassical angular momentum conservation,

m j = (R sin χ j )k j = (R + ∆ j )k out,j , (29)
where k out,j = (2m * / 2 )(E j -V 0 ). Then, the distance ∆ j becomes

∆ j = m j / (2m * / 2 )(E j -V 0 ) -R, (30)
and it is easy to see that E j = V eff (R + ∆ j ). The value of the width ∆ j for the chosen mode (m j , k j ) is confirmed numerically by evaluating the maximum of the emission Husimi function H E (D 0 = 1.5R, h) of the eigenmode, as depicted by the red solid line in Fig. 6. Let us note that the emission Husimi function of the eigenmodes does not depend on the parameter α because of their rotational symmetry.

Taking into account all these resonance modes, one can reasonably expect the width ∆ to be predicted by summing the contributions of each resonance component weighted by an appropriate factor, and thus be expressed by the average value over all the ∆ j 's, where γ j is the decay rate of the mode, appearing in the imaginary part of the energy E j = ω j -iγ j /2. Performing the sum over the 608 relevant resonant modes in the case of the wave packet B, Eq. (31) leads to ∆ ≃ 0.30 which is consistent with the previous numerical observation ∆ ≃ 0.324.

∆ = j |b j | 2 γ j ∆ j j |b j | 2 γ j , (31)

C. Time delay of tunneling wave packet

In this subsection, we shall provide a qualitative explanation of the rotation angle α T of the tunneling wave packet by investigating the time delay of 1D tunneling wave packets evolving in a 1D effective potential.

First, we start with a simple rectangular barrier as shown in Fig. 7 (a) and defined as

V rec =      0 if x < x a ; V max if x a ≤ x ≤ x b ; V min if x > x b , (32)
where we shall choose V max = 100, V min = 60, x a = 0 and x b = 1 for the following numerical computations.

We consider an incident wave packet composed of plane wave components,

Ψ(x, t) = dk W (k, k 0 )e i(k(x-xa)-E(k)t/ ) , (33)
where the window function is given by W (k, k 0 ) = e -(k-k0) 2 /σ 2 and the central energy of the incident wave packet is E 0 = ( k 0 ) 2 /2m * . This incident wave packet would hit the right side of the barrier (x = x a ) at t = 0. We can get the tunneling (transmitted) wave packet appearing beyond the barrier by inserting the transmission coefficient T (k) [START_REF] Liboff | Introductory Quantum Mechanics[END_REF],

Ψ(x, t) = dk W (k, k 0 )T (k)e i(k ′ (k)(x-x b )-E(k)t/ ) , (34)
where k ′ (k) = (2m * / 2 )(E(k) -V min ). This tunneling wave packet is set to be located on the right side of the barrier (x = x b ) at t = 0 if T (k) is independent of k (or equally E). However, T (k) depends on k or E in general and it is responsible for the time delay as discussed below.

First, let us take a narrow window function W (k, k 0 ) with a small σ to construct a wave packet Ψ H (x, t), choosing E 0 = ( k 0 ) 2 /2m * = 80 such that its central energy lies in the range 60 < E < 100 (see Fig. 7 (a)) where the barrier has a finite but nonzero width. The weighting function W (k, k 0 )T (k) of the tunneling wave packet in Eq. ( 34) is shown as a function of E (the red line) in Fig. 7 (b). Writing the transmission coefficient as T (E) = |T (E)|e iϕT and then performing the stationaryphase approximation on the integral of Eq. (34), one unravels the semiclassical dynamics of the transmitted (nearly Gaussian-like) wave packet and one extracts, similarly to Eq.( 23) for the reflected case, a time delay [START_REF] Bracher | [END_REF] δt

T = dϕ T dE . (35)
The phase ϕ T is shown to increase roughly from -0.5 to 1.5 in the energy range 60 < E < 100 as depicted by the black solid line in the lower panel of Fig. 7 (b) and a linear approximation yields the time delay δt T ≃ 2/40 = 0.05, which is greater than the time delay of the reflected wave packet, δt R ≃ 3/100 = 0.03 in this case (see the black dashed line in the lower panel of Fig. 7 (b)).

The evolution of the tunneling wave packet is shown in Fig. 7 (c). The time at which the peak of the tunneling wave packet starts to appear at x b is t = 0.045 and is consistent with the previous δt T ≃ 0.05 estimated from the variation of ϕ T . Now let us consider another wave packet Ψ L (x, t) with a lower energy (E 0 = 50) and a broader window function with a larger σ. On contrary to Ψ H , the weighting function W (k, k 0 )T (k) of this tunneling wave packet does not have the shape of a Gaussian as shown by the blue line in Fig. 7 (b). Nevertheless the time delay, evaluated from the positions of the peak of the tunneling wave packet (see Fig. 7 (d)), is similar to the previous Gaussian-like packet, i.e., δt T ≃ 0.045. The average energies of the transmitted part of Ψ H and Ψ L , evaluated from the velocity of the peak of the wave packets in Fig. 7 (c) and Fig. 7 (d), respectively, are depicted by the red and blue vertical dotted lines on Fig. 7 (b). One can note that, for both energies, the slope of the phase ϕ T (lower panel in Fig. 7 (b)) is very similar and therefore, according to Eq. ( 35), fully consistent with the comparable time delays δt T observed.

One can finally obtain the time delay of the transmitted 2D wave packet studied in Sec. V. A by considering now the effective potential from Eq. ( 4) modified as shown in Fig. 4 (a), taking m = 17 and V 0 = 100 which is nothing but a rescaling of the potential we have been using in Sec. IV B, since m 2 /V 0 = 120 2 /5000 ≃ 17 2 /100 ≃ 2.9. We take a narrow window function W (k, k 0 ) with E 0 ≃ 118 as shown in Figs. 8 (a,b). Taking into account the shift due to T (E), the right side of the barrier is then evaluated at x b ≃ 2.7. The transmitted coefficient T (E) is again determined by matching the plane wave solutions and the Hankel functions on both sides of the barrier as already explained in Sec. IV B. The time evolution of the tunneling wave packet is shown in Fig. 8 (c). Once again, the time delay can be determined by the maximum point of |Ψ(x b , t)| 2 as illustrated in Fig. 8 (d). If we assume x b = 2.7, we finally obtain the time delay of δt T ≃ 0.14. Note that this is larger than the reflected time delay δt R ≃ 3/100 = 0.03. Based on the previ- ous calculations for the rectangular barrier, we can also assume the ratio δt T /δt R ≃ 4.6 to be almost constant even when the tunneling wave packet is not Gaussianlike. Then, one can apply this ratio to our effective potential Eq. ( 4) with m = 120 and V 0 = 5000. Therefore the expected time delay of the tunneling wave packet in this case would be δt T ≃ 4.6δt R = 4.6 × (3/5000) ≃ 0.00276. This is consistent with the time delay δt * = 0.227s 0 = 0.227/90 ≃ 0.0025 deduced from the dynamics of the 2D tunneling wave packet in Sec. V A.

VI. HIGH-ENERGY WAVE PACKETS

Finally, we consider the high-energy wave packets C and D. The central energy of the wave packet C with (m 0 , k 0 ) ≃ (140, 122) and χ ≃ 35.0 degrees is chosen to correspond to the top of the barrier of the effective potential V eff defined by Eq. ( 4) (see Fig. 1 (b)) such that the incident wave packet would be roughly equally distributed into a transmitted and reflected part after the bounce as shown in Fig. 9 (a) (61.6% of incident wave packet is reflected). Like for the wave packets A and B, the reflected wave packet exhibits a GH shift and an increased reflection angle χ R which enhances as the central energy of the wave packet gets larger. Besides, the transmitted wave packet exhibits a transmission angle of 64.7 degrees, which is a considerable deviation from the tangential direction (the classical expectation of the transmission angle 90.0 degrees) (see Fig. 9 (c)). This corresponds to the Fresnel filtering effect on the transmitted wave packet, well-known in optics [25].

When the energy becomes even higher like the wave packet D with (m 0 , k 0 ) ≃ (140, 140) and χ ≃ 30.0 degrees, almost the whole wave packet (about 91.6 %) would be transmitted as shown in Fig. 9 (b) and (d). In this deep semiclassical regime, the wave packet follows the classical trajectory of a particle without any appreciable lateral shift.

VII. SUMMARY

We have studied the dynamics of wave packets in the 2D circular step potential. The time evolution of the wave packets has been obtained from the interference between the time oscillations of eigenmodes of the system, i.e., bound modes and resonance modes. At low energy, the wave packets exhibit a GH shift along the circular boundary when they are reflected by the curved step potential. We have shown that the GH shift is well explained by the variation of the phase of the reflection coefficient which is obtained using a modified 1D effective potential deriving from the 2D system. The good agreement with the numerical simulations justifies our assumption that only the shape of the barrier region is relevant to describe the time delay due to reflection in a 1D potential. We have also demonstrated that the tunneling wave packet is characterized by a finite distance ∆ from the curved boundary by investigating its dynamics with the help of the emission Husimi function. The gap ∆ can be estimated by averaging over each individual resonance contributions. We believe that this wave packet approach will furnish useful materials in order to investigate the natural connections between classical and quantum mechanics.

FIG. 1 :

 1 FIG. 1: (Color online) (a) The effective potential V eff (r) defined by Eq. (4) for m = 120 and V0 = 5000 is depicted by the blue line. The yellow and purple color correspond to the energy ranges where the bound and tunneling modes lie, respectively. (b) For the same V0, the eigenmodes are represented in the (m, Re(k)) plane. The central values of m and k of four wave packets (A-D) are marked by large solid dots, i.e. (m0, k0) = (75, 75) (red), (120, 90) (blue), (140, 122.065) (green), and (140, 140) (yellow), respectively. The tiny orange dots around B mark the relevant eigenmodes used to reconstruct the wave packet B (see Sec. III).

FIG. 2 :

 2 FIG. 2: (Color online) (a) The initial Gaussian wave packet |ΨI (x, y)| 2 with (m0, k0) = (120, 90) is plotted with the 2D potential in the (x, y)-plane. (b) Comparison between the original Gaussian wave packet defined by Eq. (18), depicted by the red line, and the reconstructed wave packet (dashed green line) obtained from Eq. (19). The horizontal and vertical profiles of |ΨI (x, y)| 2 are shown. The center of |ΨI (x, y)| 2 is located at (x0, y0).

FIG. 3 :

 3 FIG. 3: (Color online) The time evolution of the wave packets is shown in the (x, y)-plane. The wave packets A and B are plotted in (a) and (b), respectively, at t = 0.5 and 1.5. The thick vertical yellow line corresponds to the y-axis while the white curved line depicts the boundary of the circular cavity.In the panel (c) and (d), the red rectangles are the expected positions of the wave packets A and B at various times obtained from Eq. (24) with 1 % and 2.5 % increased reflection angles, respectively, and show a very good agreement with the average positions xav(t) denoted by the solid blue dots.

FIG. 4 :

 4 FIG. 4: (Color online)The effective potential V eff (r) (Eq. (4)) is modified according to the panel (a), from left to right, in order to define the phase change ϕR upon reflection. The effective potential is plotted for m = 120 and V0 = 5000. On the panel (b) the phase change ϕR (red solid line) is plotted with respect to the energy and compared to the phase φ of the 1D step potential case (see Eq. (22)) with V0 = 5000. The phases are nearly the same except for higher energy as shown in the inset.

FIG. 5 :

 5 FIG. 5: (Color online) Tunneling transmission of the wave packet B. (a) Schematic illustration of the coordinates (D, h) and the rotation angle α for the emission Husimi function HE(D, h; α). (b) The strength of tunneling transmission f (α) defined by Eq. (28) is evaluated at t = 7.5 (red dotted line). The vertical green dotted line indicates the maximum of f (α). (c) The distance ∆ of the normal tunneling wave component is plotted as a function of the angle α at two different times t = 7.5 and t = 10 (respectively the red solid dots and the open blue circles). The vertical green dotted line marks the value αT = -0.045 where the curves intersect. (d) The normal tunneling wave component is plotted at t = 7.5 and t = 10 for the particular value αT = -0.045. (e) The average positions xav of the normal tunneling wave component are shown in the (x, y)-plane. The position at t = 1 (solid red triangle) is extrapolated from the observed positions at t = 7.5 and t = 10 (solid red circles).

1 (

 1 µπ) 1/4 e -(h ′ -h) 2 /2µ with µ = (0.15) 2 /2. This emission Husimi function probes the intensity of the wave component normally incident to the projection line. Any point (x, y) is related to (D, h) at an angle α asx = D cos α -h sin α, y = D sin α + h cos α. (27)For a given value of α, we can track the tunneling wave component, normally incident to the h-axis, by plotting H E (D, h), as shown in Fig.5(a) and (d). We can then obtain the average position (D 0 , h 0 ) of the emission Husimi function in the (D, h) plane and the strength of the normal tunneling component defined asf (α) = dD dh H E (D, h; α)(28) at a given time t. The strength f (α) evaluated at t = 7.5 is shown as a function of α in Fig. 5 (b) and the plot exhibits a maximum at α T ≃ -0.045. This angle α T indicates the most probable direction of the tunneling transmission (see Fig. 5 (d)). This most probable tunneling direction can be confirmed with the plot of ∆ = h 0 -R. The values of ∆ at two distinct times, t = 7.5 and t = 10 are shown as a function of α in Fig. 5 (c) and it is clear that we can get the time-independent value ∆ ≃ 0.324 from the crossing point of the two curves which corresponds to the tunneling direction α T . The figure 5 (d) shows the evolution of the normal tunneling wave component by plotting H E (D, h; α T ) at t = 7.5 and 10. The emission Husimi function spreads out as the time increases, remaining ∆ invariant. The negative value of the angle α T implies that the tangential tunneling wave component appears beyond the barrier after some time delay. The initial position of the tangential tunneling wave component, i.e. (x T , y T ) ≃ ((R + ∆) sin |α T |, (R + ∆) cos α T ) ≃ (0.105, 2.321), is marked by the open red circle in Fig. 5 (d). The transmitted wave component propagates freely and from its average positions, (x, y) = (2.993, 2.192) at t = 7.5 and (4.146, 2.14) at t = 10, we can estimate the position (x, y) = (-0.005, 2.326) at the time t = 1.

FIG. 6 :

 6 FIG. 6: (Color online) The resonance mode (m, k) = (120, 113 -i1.57 × 10 -6 ) is shown in the upper part of the picture. In the lower part, the blue line is the corresponding effective potential with V0 = 5000 and the red line represents the emission Husimi function HE(D0, h) of the mode with D0 = 1.5R, whose maximum gives the individual gap ∆j . The horizontal dashed line superimposed on the potential denotes the energy Ej of the resonance mode.

FIG. 7 :

 7 FIG. 7: (Color online) Dynamics of tunneling wave packets in a rectangular barrier. (a) The potential shape depicted here is given by Eq. (32) with Vmax = 100, V0 = 60, xa = 0 and x b = 1. The red (blue) arrow approaching the barrier denote the energie of the incident wave packet ΨH (ΨL). The average energy of tunneling wave packets are marked by the arrows beyond the barrier. (b) The upper panel describes the weighting distribution W (E)|T (E)| of the plane wave components of the tunneling wave packets ΨH and ΨL (red and blue line, respectively) as a function of the energy. The lower panel depicts the phase ϕT of the transmission coefficient T (E) (black solid line). The phase of the reflection coefficient (black dashed line) is also plotted for comparison. The time evolution of the tunneling wave packets |ΨH| 2 and |ΨL| 2 is shown in (c) and (d), respectively.

FIG. 8 :

 8 FIG. 8: (Color online) The dynamics of the tunneling wave packet with an average energy E0 ≃ 118 (horizontal orange arrow in (a)) is investigated in the modified effective potential shown in (a) and defined in Fig. 4 for the parameters m = 17 and V0 = 100. In the panel (b), the distribution of the plane wave components of the incident wave packet is plotted (red line) with respect to the energy and the maximum of the function is denoted by the vertical orange dotted line. The time evolution of the square modulus of the wave packet is depicted in (c). The time variation of |Ψ(x, t)| 2 at x = 2.7 and 2.8 is displayed in (d) and the time range delimited by the two maxima (light red window) corresponds to the time delay of tunneling wave packet.

FIG. 9 :

 9 FIG. 9: (Color online) Time evolution of the high-energy wave packets. The wave packets C and D are respectively plotted in (a) and (b) in the (x, y)-plane at t = 0.6 and 1.4 and the average positions xav(t) at various times are denoted by red solid dots in (c) and (d) for the wave packet C and D, respectively.

VIII. ACKNOWLEDGMENTS

This research was supported by Kyungpook National University Research Fund, 2013, and Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education(No. 2013R1A1A2065357). We also acknowledge support by the Deutsche Forschungsgemeinschaft within the Forschergruppe 760 "Scattering Systems with Complex Dynamics."