Mathieu Guidere

The Verbal Psychotrauma

Keywords:

Psychotrauma is the result of life-threatening attack or a natural disaster. In many expressed in different ways, the most expression is rarely taken into account a better understanding of the mental method for analyzing traumatic expressions

The

  Journal of Applied Research in Human & Social Sciences (JARHSS) N°2 -Pages 7-20 -ISSN 2677-7193. www.journal-of-applied


Marie Masson
email: marie.masson@orange.com

Zwi Altman
email: zwi.altman@orange.com

Eitan Altman
email: eitan.altman@sophia.inria.fr

Coordinated scheduling based on automatic neighbor beam relation

Keywords: Beam relations, ANBR, Massive MIMO, Coordination, Multi-user scheduling, Interference management, 5G --

We consider a network with Massive Multiple Input Multiple Output (M-MIMO) base stations using a Grid of Beams (GoB) for data and control channels. 5G allows to establish interference relations between beams of neighboring cells. Such relations can be used to automatically generate a beam relation matrix, denoted as Automatic Neighbor Beam Relation (ANBR) matrix that can be very useful for optimizing different resource allocation processes. This paper shows how the ANBR matrix can be used to coordinate scheduling of neighboring cells with a small amount of information exchange. The coordination is performed by judicially muting or reducing the bandwidth of certain beams in the process of Multi-User (MU) Proportional Fair (PF) scheduling. Numerical results show how the coordination approach can bring about significant performance gain.

I. INTRODUCTION

M-MIMO is among the pillars of 5G technology that allows to significantly improve user rates, system capacity and Energy Efficiency (EE) [START_REF] Björnson | Massive MIMO networks: Spectral, energy, and hardware efficiency[END_REF]. The concept of GoB has been introduced to allow beamforming of control channels that are used among others to transmit synchronization signals and to broadcast system information to allow initial access and mobility procedures [START_REF] Vook | MIMO and beamforming solutions for 5G technology[END_REF]. Beam sweeping is used in conjunction with GoB by switched rapidly the beams one by one in a manner to cover the entire cell surface. Similarly, data channels can be transmitted using the beams of the GoB. In this case the GoB can be seen as a predefined codebook of beams that can be selected by a MU scheduler to serve users.

Beam management is an important feature introduced in 3rd Generation Partnership Project (3GPP) for 5G networks. It is part of the New Radio (NR) Automatic Neighbor Relation (ANR) and is considered as a central Self-Organizing Network (SON) function. NR ANR allows to automatically establish different types of relations involving gNodeBs (gNBs) and/or beams [START_REF] Ramachandra | Automatic neighbor relations (ANR) in 3GPP NR[END_REF]: (i) gNB to gNB relations that consist of establishing connectivity over the Xn interface between neighboring gNBs. Such relations are necessary to support mobility, load and traffic sharing or multi-connectivity and were already standardized for 4G networks [START_REF] Dahlen | Evaluations of lte automatic neighbor relations[END_REF]; (ii) gNBs to Beam relations; and (iii) beam to beam relations, intraand inter-cell.

The feature of automatically establishing relations between beams is denoted as ANBR. The association of user traffic and Quality of Service (QoS) information with the ANBR feedback in relation with the GoB enable the exploitation of the beam level spatial resolution to further optimize resource management functions.

This paper investigates the use of ANBR to optimize MU-Multiple Input Multiple Output (MIMO) scheduling. We assume that the ANBR feature is available and provides a binary static matrix, denoted as ANBR matrix, with non-zero elements representing beam relations. The way to define relations between beams is not standardized. It can be derived as in classical ANR, e.g. by calculating the average interference users of beam i experience from beam j of a neighboring cell during a long period of time, and comparing it to a predefined threshold. We show how the ANBR matrix can be used to coordinate the MU schedulers of neighboring cells in a manner to minimize collisions between potentially interfering beams, while taking into account the dynamics of the traffic.

Traditional approaches for collaboration between adjacent cells with M-MIMO deployment such as Coordinated Multipoint transmission (CoMP) [START_REF] Irmer | Coordinated multipoint: Concepts, performance, and field trial results[END_REF] require signaling and computation in both the PHY and the MAC layer and are thus demanding an intensive signal processing and the exchange of important amount of information. They also have strict requirements on the back-haul capacity. In the proposed approach, the coordination is performed at the MAC level and requires little information echange and processing power.

The paper is organized as follows. Section II presents the system model. Section III develops a coordinated MU scheduling solution relying on the ANBR. The architecture supporting the solution is also described. Numerical results are presented in Section IV followed by concluding remarks in Section V.

II. SYSTEM MODEL

Consider a hexagonal network with tri-sectoral sites, each with M-MIMO system in the Down Link (DL). GoB with beam sweeping is used for initial access (to attach users to the Base Station (BS)) and synchronization. Data transmissions use the (fixed) beams of the GoB, and the choice of users being served is made by the MUscheduler.

Each sector (BS or cell) m is equipped with a M-MIMO antenna with N = N x × N z radiating elements, and is serving N m mobiles, each with one receiving antenna (the beam generation and antenna modeling is explained in details in [START_REF] Tall | Multilevel beamforming for high data rate communication in 5G networks[END_REF]). Figure 1 presents the coverage area provided by GoBs of two adjacent cells for the case of N x = N z = 16. The color code is used to simplify visualization and has no physical meaning. Fading is removed for clarity. It is noted that the zeros of the beam radiation patterns in the azimuth axis close to the BSs is due to the zeros in the beam radiation patterns and the fact that reflections and fading are not omitted. It is recalled that beams for control channels are activated one by one (via beam sweeping) whereas several beams for data channels can be scheduled simultaneously.

We consider two BSs (i.e. macro cells) m and m which interfere each other, each of which having a GoB denoted by B m and B m respectively as shown in Figure 1. We assume that a cell can serve up to K users in a time slot, with at most one user per beam b. The users are selected according to a PF criterion.

Assume that the BS serves k users at a given time slot with k ≤ K and denote by p max the maximum transmission power of the BS. The power p m that the BS m transmits to user u, equals

C m,b u = p m |h m b (u)| 2 . ( 1
)
|h m bu (u)| 2 is modeled as the product between the pathloss, the antenna gain G m bu (u) of the beam b u serving Fig. 1: GoBs of two neighboring cells projected on the surface user u and measured at the direction of user u, and the fast fading term Z(u). The latter is modeled as a realization (per user) of a Nakagami distribution, and can be parameterized for different propagation environment [START_REF] Charash | Reception through Nakagami fading multipath channels with random delays[END_REF].

|h m bu (u)| 2 = c d γ u,m G m bu (u)Z(u) (2)
where c and γ are constants that depend on the type of environment. The interference generated by a cell m on u is written as the sum of interferences from its active beams I b u :

I m u = b ∈B m I b u (3)
where

I b u = p m u ∈b ,u =u |h m b u (u)| 2 (4)
The Signal to Interference plus Noise Ratio (SINR) of a user u attached to cell m is written as:

S m u = p m |h m bu (u)| 2 m I m u + σ 2 (5)
A full buffer traffic model is assumed. The baseline scheduler is based on PF without coordination.

III. COORDINATED SELF-ORGANIZING SCHEDULING

A. Architecture framework

The architecture supporting the ANBR based coordination is shown in (Figure 2). A centralized ANBR SON function is deployed at the management and orchestration plane as an application (but can be implemented in a distributed manner as well). It provides a static matrix A m,m for any two cells m and m , with elements A m,m b,b . For simplicity of notations we omit the superscripts m and m from the matrix A in the rest of the sequel. A b,b = 1 if the beam b of cell m and beam b of cell m interfere each other (as mentioned in Section I), and 0 otherwise. The matrix A is calculated and updated (not often) according to the operator policy and is considered here as constant.

The Distributed Unit (DU) (i.e. Radio Link Control (RLC), Media Access Control (MAC) and Physical (PHY)-High control protocol stack), hosts a new functional block denoted in Figure 2 as D-ANBR. The D-ANBR receives and processes measurements from the Radio Units (RUs). It dynamically updates the beam relations' matrix A and keeps only those relations corresponding to the present traffic distribution (as described in the next subsection). The resulting sparser matrix, Q, is transmitted to the schedulers of BSs m and m and is used to coordinate them. The time scale for generating the matrix Q is that of the traffic dynamics, e.g. of the order of a second.

It is noted that the introduction of new control, Radio Resource Management (RRM) or machine learning algorithms is currently studied in different standardization fora. As an example, the ORAN Alliance is defining new interfaces and protocols to support deployment of such algorithms in different network nodes [START_REF] Alliance | O-RAN:towards an open and smart ran[END_REF].

U serv b,m = {u ∈ m|b = arg max b∈Bm C m,b u } (6)
We define two indicators used to generate the matrix Q using A. The first, A 1 (b, u ), indicates whether user u achieves low Signal to Interference Ratio (SIR) that is below a predefined threshold γ th and is thus likely to experience strong interference from beam b:

A 1 (b, u ) = 1 if C m u I b u < γ th 0 otherwise (7)
For clarity of notation the beam b serving user u is not included in [START_REF] Charash | Reception through Nakagami fading multipath channels with random delays[END_REF]. Denote by U int b the set of users u ∈ m for which A 1 (b, u ) = 1, namely the set of users that could benefit the most from reducing the interference from beam b:

U int b = {u ∈ m |A 1 (b, u ) = 1}. (8)
The cardinality of U serv The second indicator, A 2 (b), is used to verify whether the ratio between the number of users that beam b interferes and the number of users it serves is above a threshold η th . The rational is that the benefit from muting or limiting the allocated resources to a beam during a Transmission Time Interval (TTI) increases with the amount of users it interferes and decreases with the amount of users it serves.

A 2 (b) = 1 if n int b n serv b > η th 0 otherwise. ( 11
)
The matrix element Q b,b of Q is defined as follows:

Q b,b = A b,b × 1 {A2(b)=1} × 1 {A2(b )=1} (12)
The threshhold values of η th and γ th are determined using a simple optimization procedure (see details is Section IV).

The rationale for (12) is the following: consider the case where coordination is based on constraining the MU schedulers, namely not to schedule users served by beams b and b in the same TTI for which Q b,b = 1 (see Algorithm 2). From eq. ( 12), the coordination is performed if both users served by b and b will benefit from coordination, namely local QoS loss is avoided even if globally performance is improved. The algorithm for generating the matrix Q is given by Algorithm 1.

Algorithm 1 Code Generation

Input: A Init: Q =[0], matrix of the same dimension as A for all couples (b, b ) for which A b,b = 1 do n int b ← |U int b |, n serv b ← |U serv b | n int b ← |U int b |, n serv b ← |U serv b | Q b,b = A b,b × 1 {A2(b)} × 1 {A2(b )} end for
The MU scheduling uses a known beam selection feature known as beam skipping technique that is applied independently in each cell. It allows to reduce intra-and inter-beam interference and to improve users' rates. It is noted that the ANBR based coordinated scheduling is independent of the beam skipping feature and can be applied without it. Each user is attached to the beam of the GoB achieving the best SINR. The user attachment provides certain spatial information that is exploited by the scheduler to mitigate interference: (i) by avoiding scheduling two users attached to the same control beam and (ii) by avoiding scheduling two users attached to adjacent beams of the GoB, both in the same TTI.

We first present a time based ANBR coordination scheme, denoted for sake of brevity as time-ANBR scheme. In this scheme, coordination is achieved by muting of certain beams during certain TTIs as explained below. For sake of simplicity, delay has been ignored but can be easily incorporated into the coordiation scheme.

Denote by R b u,t M +1 the instantaneous rate of a user u ∈ m when it is scheduled at t M +1 , by R b u,t M -its average rate computed over a time window [t M -T , t M ], and by d > 0 a number chosen as small as possible. Denote by U candidates the set of users that can still be scheduled and by U K -the set of users already selected for scheduling, both at t M +1 .

Consider next the scheduling algorithm of cell m (or m ) (see Algorithm 2). In the initialization phase, U candidates contains all the users attached to m (or m ). The scheduler ranks the users in U candidates with respect to a PF criterion, namely

R b u,t M +1 R b u,t M +d .
All beams b ∈ B m for which Q b,b = 1 are muted at an even TTI, whereas the beams b ∈ B m are muted at an odd TTI. The users of a muted beam are removed from U candidates . The scheduler selects the top-ranked candidate. Then, it removes the selected users attached to the adjacent beams from the candidate list (following the beam skipping scheme). We repeat the above two operations until K candidates are selected or until the set of candidates is empty.

Algorithm 2 ANBR-assisted MU Scheduler

Input : Q Init: U K = {}, U candidates = {u ∈ m} if TTI is even then for all b for which Q b,b = 1 do Remove u ∈ b from U candidates end for end if while |U k | < k or U candidates = ∅ do u select ← argmax u∈U candidates R b u,t M +1 R b u,t M +d U k ← U k ∪ u select
Remove users attached to the beam of u select and to the adjacent beams from U candidates end while

The second coordination scheme is denoted as frequency based ANBR coordination scheme and is denoted for sake of brevity as frequency-ANBR scheme. This scheme is similar to the time-ANBR namely instead of consecutively muting beams for which Q b,b = 1, we allocate to these beams half of the available nonoverlapping resources. It is noted that in spite of being very simple, both time-and frequency-ANBR coordination schemes achieve high performance with little computational efforts.

Lastly, the generalization of the coordinated scheduling to the case where beams from three cells interfere with each other is straightforwards. Sparse three dimensional matrices A and Q need to be generated, and nonoverlapping resources (e.g. Physical Resource Blocks (PRBs) in the frequency-ANBR) should be allocated to the beams. One should bare in mind that significant co-located interference from several cells should be minimized in the cell-planning phase and not by the scheduler.

IV. NUMERICAL RESULTS

A. Simulation scenario

Consider a network comprising 19 sectors (cells): a central sector and two tiers of 18 neighboring sectors (6 sites located on a hexagonal grid) surrounding it. Each base station is equipped with a M-MIMO antenna as described in Section II. The central sector and one of its direct neighbors denoted hereafter as cell 1 and cell 2 (respectively upper left and lower right in Figure 1), implement the coordinated scheduling. The reference scenario with no coordination serves as a baseline. It implements a PF based MU-scheduler with the beam skipping feature. Both the time-and frequency-based ANBR coordinated scheduling are simulated and compared to the baseline scenario. The simulation parameters are summarized in Table I. The traffic distribution of cells 1 and 2 is shown in Figure 4. A red rectangle surrounds a hotspot zone with high traffic located around the cell edge area of the two cells. Each cell has 35 users in the hotspot area and 10 users in the rest of the cell, drawn according to a uniform distribution in each zone. The users' color code in Figure 4 is the following: red and blue squares are the users outside of the hotspot and belonging to the cell 1 and 2 respectively, green and yellow are the users of cell 1 and 2 located in the hotspot.

TABLE I: Network and Traffic characteristics

The best values for the thresholds η th and γ th are determined by means of an exhaustive search. We define a uniform grid of 10 × 10 points (η th ,γ th ), with η th varying from 1/10 to 1 and γ th -from 1 to 10. For each point of the grid we compute the Mean User Throughput (MUT) gain using the time-ANBR coordination scheme as depicted in Figure 3. The gain increases with the decrease in η th while the MUT gain is not sensitive to variations in γ th except for small values. Too small value of γ th results in too few mobiles that can benefit from the coordination between the two cells. Similarly, a small value of η th allows more beam relations to be included in the matrix Q and more users will participate in the coordinated scheduling. In the rest of the paper, the thresholds' values of γ th = 2.5 and η th = 1/6 are set. A gain of 105% in MUT with respect to the baseline is achieved on a plateau of 27 points, indicating little sensitivity of the thresholds (see Figure 3).

Figure 5 presents the distribution of the served -and 6 and the time evolution of the MUT for both cells is shown in Figure 7. The improvement brought about by the coordination schemes is very significant, of the order of 100%. One can see that the time-ANBR performs a Fig. 5: Number of served and interfered users per beam bit better than the frequency-ANBR based coordination. This is explained by the fact that when a beam is muted, time resources will be used by users served by other beams. In the case of frequency-ANBR, the coordinated beams are allocated non-overlapping frequency resources and hence not all the available resource are used.

In the following results, we consider the frequency-ANBR coordination scheme compared to the baseline. We divide the users into two groups as a function of their locations, namely outside and inside the hotspot area. The users' throughput are presented in the form of horizontal bars, in an increasing order of throughput values in the baseline case (in blue). The same order is kept for the coordinated scheduling case (in red) to ease comparison.

Figures 8 and 10 present the throughputs of the users outside the hotspot zone of cells 1 and 2 respectively. In cell 1 (8) certain users see their throughput grows significantly since they benefit from the cells coordina-

  pmax k . Denote by C m,b u the useful signal received power of user u from beam b of BS m, by d u,m -the distance between m and u, and by σ -the thermal noise. C m,b u can be written as a function of the channel gain h m bu (u):

Fig. 2 :

 2 Fig. 2: System architecture

  of user in the hotspot zone 35 of each cell Number of user outside the hotspot 10 zone in each cell Traffic distribution in hotspot zone Uniform Traffic distribution outside hotspot zone Uniform Service Type Full buffer, data

Fig. 3 :

 3 Fig. 3: Comparison of MUT as a function of the thresholds γ th and η th

Fig. 4 :

 4 Fig. 4: Traffic map

Fig. 6 :Fig. 7 :

 67 Fig. 6: MUT for frequency(in red)-and time(in yellow)-ANBR and baseline(in blue)

Fig. 8 : 1 Fig. 9 : 2 Fig. 11 :

 819211 Fig. 8: Throughputs of users outside the hotspot zone in cell 1


Bhaawan Gupta

Benjamin Ducharne

Gaël Sebald

Tetsuya Uchimoto

Gael Sebald

Magnetic non-destructive materials characterization of non-ferromagnetic materials using magnetic susceptibility

Keywords: non-destructive testing, magnetic force, magnetic susceptibility, paramagnetism, diamagnetism

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Micromagnetic materials characterization is a well-known field of non-destructive evaluation and testing (NDE/NDT). It is based on the interaction of magnetic domain walls (Bloch walls) with microstructure, which is similar to the interaction of dislocations with microstructure [1 -3]. Every material interacts with a magnetic field to some extent, depending on its internal composition and structure. The interactions which define the material's behavior in a magnetic field take place on the micro-and nanometer scale and are partially quantum mechanical in nature. The magnetic susceptibility  describes how a material interacts with an external magnetic field. It has been verified before that a correlation exists between micro-or electromagnetic and mechanical properties in para-and diamagnetic materials as well [START_REF] Gupta | [END_REF]. Previously, a magnetic force sensor developed [START_REF] Gupta | [END_REF][5] has shown promising potential in characterizing the materials like composites, plastics and many other non-ferromagnetic materials. Since, luxury car makers and aerospace industry are using more and more composites in manufacturing, it becomes very important to have a method to control the quality check, e.g., adhesion between two composites. However, one major drawback of this approach [START_REF] Gupta | [END_REF] was the use of cantilever which is not very useful when it comes to industrial test conditions which are more prone to vibrations. Hence, a new modification is proposed in this contribution in which the experimental set-up is changed entirely and a diaphragm is used instead of a cantilever. This system shows more stability and is less prone to vibrations. As an added advantage, the cost of the entire set-up is also reduced.

Experimental set-up and results

Previously, a capacitive distance sensor was used in order to determine the deflection of a cantilever carrying a sensing magnet [START_REF] Gupta | [END_REF]. In the new proposed method, the authors use a plastic film thick enough to hold the sensing magnet with one layer of reflector on top for the laser sensor to sense the displacement (Figure 1). The sheet used to hold the sensing magnet is PVDF in order to be strong enough to hold the magnet and at the same time thin enough to sense some stretching due to the force on the magnet when exposed to the material under test. The thickness of the diaphragm used in this set-up is 500 microns. The laser sensor used in this method is an MTI instruments device whose resolution is 0.5 μm/mV. The output voltage change of the laser sensor electronics when the sensing magnet is approaching a material is proportional to the displacement of the magnet. When investigating practical materials such as plastics or aluminum, the deflections are so small that the distance between sample and magnet changes to very small extent. In this case, the deflection (or voltage change) is almost proportional to the susceptibility of the material. The output of the laser sensor is calibrated in terms of force. Figure 2 depicts the initial results with aluminium and graphite sample over varied distances. Many other non-ferromagnetic materials characterization will be presented in the extended version of this article to verify the efficiency of the sensor.

Conclusions

In this study a new experimental set-up is proposed for the characterization of the nonferromagnetic materials non-destructively. A diaphragm is used to avoid any vibrations which are a major possibility in the industrial testing conditions. Additionally, the cost of the sensor is also reduced to make it easily available for the industry. In the extended version of the article, many different non-ferromagnetic materials will be characterized and imaged as per their magnetic susceptibility.

Acknowledgements

The authors are thankful to JSPS Core-to-Core Program, A. Advanced Research Networks for `International research core on smart layered materials and structures for energy saving' supporting this research. The Institute of Fluid Science from Tohoku University is also gratefully acknowledged for providing support through the General Collaborative Program Figure 1: Experimental set-up schematic Figure 2: Initial results project # J17I065. The authors also acknowledge the support of LGEF at INSA Lyon for providing the facility to conduct these experiments.


Michèle Colomb
email: michele.colomb@cerema.fr

Pierre Duthon
email: pierre.duthon@cerema.fr

Frédéric Bernardin
email: fredredic.bernardin@cerema.fr

Spectral reflectance

teaching and research institutions in France or abroad, or from public or private research centers.

Abstract-A wide range of sensors are developed for autonomous vehicles in order to ensure safe driving conditions by detecting any relevant objects in the road environment. The sensors, cameras, Lidars and radars, cover a wide range of wavelength, ranging from visible, Near InfraRed, Short Wave InfraRed and beyond for high frequency radars. For Lidars technology, a single wavelength scans the environment. To get the most efficient detection, the choice of the wavelength must be optimized. Knowledge of the spectral reflectance of road objects at the chosen wavelength is a key point. An original application of reflectance measurement technics is proposed to characterize road objects in the range of wavelength between 350 to 2450 nm. This investigation could lead to the definition of reference tests including a physical characterization of the environment.

I. INTRODUCTION

Intelligent transportation systems are facing an increase of sensors development to cope with the challenge of driving automation. The future intelligent systems of the autonomous vehicles must provide two major functions: navigation and obstacle detection in all possible use cases, ordinary but also in the most critical environment. The challenge is to provide these functions with an increased efficiency compared to the current situation where humans are at the controls of vehicles. Indeed, there is a wide technological offer that is proposed as cameras, Lidars or radars, in order to take the relevant visual informations to ensure the driving task instead of the driver. Until now, the driver has been monitoring the environment and taking the relevant information for the driving task in the visible range. With the development of autonomous driving, for automation level 3 to 5 [START_REF]J3016 -Taxonomy and Definitions for Terms Related to On-Road Motor Vehicle Automated Driving Systems[END_REF] the system will have to monitor the environment, and then to guarantee the driving safety [START_REF] Shladover | Technical Challenges for Fully Automated Driving Systems[END_REF].

Cameras are operating in the visible range, Near InfraRed, NIR, Short Wave InfraRed, SWIR or Long Wave InfraRed, LWIR. Relevant informations are taken from the image of the road scene, that is obtained by the reflection or emission of all objects in the considered wavelengths. It is then required to investigate what is the effect of the range of wavelength on the detection of road signaling equipment as well as on the detection of relevant fixed obstacles or of vulnerable road users. The very auguring Lidars technology is using a single wavelength to scan the environment. What is the best wavelength to use to get the most efficient sensor? Available sensors are operating at 905 nm, some studies discuss the potential benefits of using new systems operating at 1550 nm when used in foggy conditions [START_REF] Kutila | Automotive LiDAR performance verification in fog and rain[END_REF]. When combining technologies, it is also useful to know the reflectance properties of all the objects in each considered wavelength. As a single sensor is not be able to detect all objects in all conditions, fusion methods for Lidars-camera setups are proposed to optimize the object detection by computer vision [START_REF] Bijelic | Seeing Through Fog Without Seeing Fog: Deep Sensor Fusion in the Absence of Labeled Training Data[END_REF].

Multispectral characterization of road objects is made possible thanks to a spectroradiometer. Reflectance spectroscopy in the visible and NIR range has been used for years as a nondestructive tool for evaluating soil properties and vegetation studies [START_REF] Viscarra Rossel | Visible, Near Infrared, Mid Infrared or Combined Diffuse Reflectance Spectroscopy for Simultaneous Assessment of Various Soil Properties[END_REF], [START_REF] Ghulam | Characterizing Crop Responses to Background Ozone in Open-Air Agricultural Field by Using Reflectance Spectroscopy[END_REF]. Using the same physical concept, an original application of reflectance measurement technics is implemented in the Cerema Adverse Weather platform [START_REF] Colomb | An innovative artificial fog production device improved in the European project FOG[END_REF], [START_REF] Duthon | Methodology used to evaluate computer vision algorithms in adverse weather conditions[END_REF]. This paper will focus on the measurement of spectral reflection of different road objects in the range from 350 to 2450 nm. The choice of road targets is made according to the critical use cases studied for the development of autonomous vehicles and for their safety. In this first investigation, the measurements concern signaling equipment, fixed obstacles, or vulnerable road users, as pedestrians, who must be detected at a long distance, to guaranty the safety of autonomous driving.

II. STUDY OF CRITICAL USE CASES AT CEREMA ADVERSE WEATHER PLATFORM A. Critical use cases

In order to evaluate the performance of autonomous vehicles, various use cases must be considered. The identification of relevant objects in the visual scene such as lane marking, or road signs is a key point for guidance and control of trajectory of the vehicle. The detection of vulnerable road users or critical objects is essential to avoid accidents and therefore demonstrate the safety of the system and allow traffic clearance on open lanes. Indeed, perception sensors as Lidars, are the entry points for this crucial information, that is transmitted to the vehicle control system. This control system will then decide of the action: e.g. emergency braking and/or avoidance manoeuvre. Thus, use cases may include various objects as:

• in urban areas, vulnerable road users, pedestrians, cyclists and cars • on motorways, lanes marking, road signs ... • on any roads, fallen objects on the road as tires, exhaust pipes...

In facts, in a natural driving situation, environmental conditions are random, and it is not possible to repeat on request the same use cases. Therefore, testing facilities are needed to repeat the same scenarios multiple times, providing guarantees for metrological control when testing sensors. Examples of use cases investigated in the Cerema Adverse weather platform are displayed in Fig. 1, for pedestrian crossing with vulnerable road users and for work zone with signaling equipments.

B. The Cerema Adverse Weather Platform

The Cerema Adverse Weather Platform, was developed to investigate all transports systems that could be affected by adverse conditions, including fog and rain [START_REF] Colomb | An innovative artificial fog production device improved in the European project FOG[END_REF], [START_REF] Duthon | Methodology used to evaluate computer vision algorithms in adverse weather conditions[END_REF]. It allows to reproduce various scenarios, as detection of vulnerable roar users or fixed obstacles, in clear conditions, night conditions, and with various ranges of fog and rain precipitations on a total length of 30 meters (Fig. 2). Dedicated to research and development, It is also open to private companies looking for a testing facility with controlled conditions. It has been used for years in partnership or collaborative projects in order to investigate various scientific topics, as humans perception in adverse conditions [START_REF] Cavallo | Distance Perception of Vehicle Rear Lights in Fog[END_REF], [START_REF] Quétard | Combined effects of expectations and visual uncertainty upon detection and identification of a target in the fog[END_REF], vision systems capabilities in fog or rain conditions [START_REF] Bernardin | Measuring the effect of the rainfall on the windshield in terms of visual performance[END_REF], [START_REF] Marchetti | Retrieving visibility distance in fog combining infrared thermography, Principal Components Analysis and Partial Least-Square regression[END_REF], [START_REF] Pinchon | All-weather vision for automotive safety : which spectral band ?[END_REF] or computer vision algorithms for objects detection [START_REF] Dahmane | The Cerema pedestrian database : A specific database in adverse weather conditions to evaluate computer vision pedestrian detectors[END_REF], [START_REF] Bijelic | Seeing Through Fog Without Seeing Fog: Deep Sensor Fusion in the Absence of Labeled Training Data[END_REF]. The physical characteristics of rain and fog produced in the platform are described in the recent study on Lidars performances in fog and rain [START_REF] Kutila | Automotive LiDAR performance verification in fog and rain[END_REF].

This platform has got a high-level instrumentation to evaluate the performance of perceptual sensors for autonomous vehicles in adverse conditions. Some weather instruments are dedicated to characterizing the atmosphere in fog or rain conditions:

• transmissiometer, for meteorological visibility in fog from 5 to 1000 m, • optical granulometer, for fog droplet size distribution from 0.4 to 40 µm, • rain gauge and a spectro-pluviometer for rainfall rate from 0.001 to 1200 mm/h. In order to characterize the road environment itself and the object properties to be perceived by humans or by artificial vision sensors as cameras or Lidars, other devices are available:

• a video-photocolorimeter for Luminance from 0.003 to 50000 cd/m 2 in the visible range, • reference cameras in visible, NIR, SWIR and LWIR,

• and a spectroradiometer for spectral reflectance of surfaces in the range of 350 to 2450 nm.

III. SPECTRAL REFLECTANCE MEASUREMENT

A. The PSR+3500 Spectroradiometer

The spectroradiometer used in the Cerema Adverse weather platform is the PSR+3500 from Spectral Evolution (Fig. 3). This type of device measures the radiance [START_REF] Howell | Thermal radiation heat transfer[END_REF]. This radiometric quantity is expressed in W/m 2 /nm/sr. The measuring range is between 350 nm and 2450 nm. The Spectral Evolution PSR+3500 spectroradiometer and its 1 ˚aperture optics.

Light energy enters the spectroradiometer through the input optics. The light entering the sensor is collimated by a mirror. Each part of the spectrum must be directed to the photodiodes detectors corresponding to its spectral range. The sensor is composed of three photodiode arrays:

• A silicium detector measures from 350 nm to 1000 nm with a resolution of 3 nm, and uncertainty of 1 to 4%. • A cooled InGaAS sensor operates in the spectral range 1000 nm to 1900 nm with a resolution of 8 nm, and uncertainty smaller than 1%. • A second cooled InGaAS sensor measures the light power of radiation between 1900 nm and 2500 nm with a resolution of 6 nm, and uncertainty of about 1%.

B. Experimental protocol

The experimental protocol consists in carrying out a relative measurement of the reflectance of the target compared to a calibrated white target of uniform reflectance. In order to measure the spectral reflectance of object it is necessary to illuminate with a calibrated light source emitting in the corresponding wavelength range, a 400 W halogen lamp is used. During these tests, this lamp is at a fixed position, and the spectroradiometer changes of position, see Fig. 4. The target is placed at a fixed distance d, at 2 meters from the halogen source. The beam emitted by the source is orthogonal to the target surface. The spectroradiometer, target and halogen lamp are in the same horizontal plan, at 1.20 meters above the ground. Measurements are made in a dark room. The spectroradiometer is first placed to obtain a measuring angle β of 10 ˚. β is the angle formed between the source and the spectroradiometer in the horizontal plane. A reference measurement is made with a white calibrated target. Then the object to be characterized is installed in place of the calibrated target. To get the angular reflectance, the spectroradiometer is placed at angle β from 10 ˚to 80 ˚in 10 ˚steps. The same protocol is repeated for the different targets.

C. Reflectance calculation

The spectral response represents the reflectance R as a function of the wavelength λ. The reflectance R of an object corresponds to the ratio of the reflected energy r to the incident energy E received by the object under study [START_REF] Howell | Thermal radiation heat transfer[END_REF]:

R = r E (1)
In the case of the measurement set up, the reflected energy r corresponds to the measurement of the radiance performed on the object at the different measurement angles. In order to calculate the incident energy, E, it is necessary, to measure a reference radiance Rad ref , on a calibrated white target. The incident energy E calculated as follows using the spectral response of this reference target R ref :

E = Rad ref R ref (2)
Thus, by combining the two equations, it is possible to obtain the spectral response of the studied object:

R = r Rad ref R ref (3)
In order to evaluate the response of the targets as a function of the measurement angle, this formula is applied independently for each measurement angle. Three reference targets have been chosen for having a uniform reflection on the range of the studied wavelength: 350 to 2450 nm. The reference calibrated targets are Zenith Polymer from Pro-Lite Technology (Fig. 5).

IV. RESULTS

A. Validation of the spectral response of calibrated targets

In order to validate the experimental protocol and reflectance calculation, we have measured the reflectance of two other calibrated targets having uniform reflectance. The spectral response measured is compared to the reference value given in the calibration certificate. Fig. 6.a and 7.a show the results of the measurements of R, the spectral reflectance, respectively for the grey calibrated target (50% reflectance) and the black calibrated target (5% reflectance).

B. Results of the measurements on road objects

Based on the use cases described before, many measurements were made on a wide variety of road objects. We will present here a selection of objects according to their reflection characteristics: retroreflective targets, almost Lambertian objects and other objects categories. The selected retroreflective objects are:

• work zone signaling cone, with grey reflective strip;

• road reflector. The selected quasi-Lambertian reflective objects are:

• asphalt pavement;

• for pedestrian clothing : dark blue polar sweater. The other selected object on road surface are:

• tire; • exhaust pipe. In order to avoid increasing the complexity of the following figures, the uncertainty is not reported, it is of the same magnitude as above. from the left to the right, the object, its spectral response and the angular distribution of the reflectance. The analysis of each object reflectance gives the following results.

In Fig. 8, the grey part of the cone is highly reflective for an angle of 10 ˚. The reflectance does not change with the wavelength. Fig. 9 shows that the reflector reflectance is decreasing with wavelength. The most reflected radiation is that of the visible. Beyond 1550 nm, it reflects only 50% of the radiation, while its spectral response is negligible at 2400 nm. It is therefore necessary to pay attention to this type of equipment, especially with Lidars using a 1550 nm laser. In Fig. 10, the asphalt is not very reflective and it has a "quasi-Lambertian" behavior. Its spectral response does not exceed 10% in the infrared, with a maximum of 7.5% in the visible range. For the dark blue polar sweater, infrared wavelengths between 800 nm and 1550 nm are the more reflecting (Fig. 11). The tire in Fig. 12 has a higher reflectance between 750 and 1800 nm and about 2000 and 2400 nm. A narrow band of lower reflectance is observed at about 1900 nm. In Fig. 13, spectral response of the exhaust pipe is increasing with the wavelengths. Indeed, the spectral response is low, 10% for the wavelength 400 nm while the spectral response for the same angle is 70% for the wavelength 2400 nm. Its shiny appearance confirms the reflectance peak at 10 ˚.

C. Discussion

As a general result, this inventory of the reflective properties of various road objects is very useful to consider, when interpreting camera or Lidars measurements dedicated to object detection. It can also be used for simulations. Some objects may be "visible", in the sense of detectable, with a high reflectance in some wavelengths, and maybe "not very perceptible", with low reflectance, in other wavelengths, e.g. the reflector on Fig. 9. These results are in favor of the development of a sensors suite with various wavelength. Then data fusion of all measurements will optimize the system and even produce redundancy for the security of the system. These results of measurements show a zone of most important reflectance between 750 and 1550 nm that is much in favor for sensors development.

V. CONCLUSION

In the context of the development of the optical sensors for autonomous driving, the definition of a reference test to evaluate the performance of the sensors is necessary. As these sensors, as cameras or Lidars, cover wavelength ranges beyond the visible, the characterization of reflection properties of objects for all wavelengths ranges, is a prerequisite for a good understanding of the factors influencing the performance i.e. the ability to detect objects. This study highlights the need to characterize the spectral reflectance of road objects at all wavelengths. This first characterization in a "clear" atmosphere could be combined with the atmospheric attenuation effect in the case of adverse situations such as rain or fog, for a complete overview of the most critical situations for driving autonomous vehicles. These original data could also contribute to computer vision modelling and simulation of sensors detection for various objects.

  Spectral reflectance characterization of road environment to optimize the choice of autonomous vehicle sensors* Michèle Colomb 1 , Pierre Duthon 1 and Frédéric Bernardin 1

Fig. 1 .

 1 Fig. 1. Examples of use cases studied in the Cerema Adverse Weather Platform.

Fig. 2 .

 2 Fig. 2. The Cerema Adverse Weather platform [9].

Fig. 3 .

 3 Fig. 3.

Fig. 4 .

 4 Fig. 4. Diagram of the experimental set-up, source, calibrated target, spectroradiometer.

Fig. 5 .

 5 Fig. 5. The reference targets grey (50%), black (5%) and white (90%) Zenith Polymer from Pro-Lite Technology.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig.[START_REF] Ghulam | Characterizing Crop Responses to Background Ozone in Open-Air Agricultural Field by Using Reflectance Spectroscopy[END_REF].Spectral response of the grey target (50% reflectance), with uncertainty.

Fig. 9 .

 9 Fig. 9. The road reflector, spectral response and angular distribution.

Fig. 10 .

 10 Fig. 10. Asphalt pavement, spectral response and angular distribution.

Fig. 11 .

 11 Fig. 11. The dark blue polar sweater, spectral response and angular distribution.

Fig. 12 .

 12 Fig. 12. The tire, spectral response and angular distribution.

Fig. 13 .

 13 Fig.[START_REF] Marchetti | Retrieving visibility distance in fog combining infrared thermography, Principal Components Analysis and Partial Least-Square regression[END_REF]. The exhaust pipe, spectral response and angular distribution.

Figures 8 to 13

 13 Figures 8 to[START_REF] Marchetti | Retrieving visibility distance in fog combining infrared thermography, Principal Components Analysis and Partial Least-Square regression[END_REF] show from the left to the right, the object, its spectral response and the angular distribution of the reflectance. The analysis of each object reflectance gives the following results.In Fig.8, the grey part of the cone is highly reflective for an angle of 10 ˚. The reflectance does not change with the wavelength. Fig.9shows that the reflector reflectance is decreasing with wavelength. The most reflected radiation is that of the visible. Beyond 1550 nm, it reflects only 50% of the radiation, while its spectral response is negligible at 2400 nm. It is therefore necessary to pay attention to this type of equipment, especially with Lidars using a 1550 nm laser. In Fig.10, the asphalt is not very reflective and it has a "quasi-Lambertian" behavior. Its spectral response does not exceed 10% in the infrared, with a maximum of 7.5% in the visible range. For the dark blue polar sweater, infrared wavelengths between 800 nm and 1550 nm are the more reflecting (Fig.11). The tire in Fig.12has a higher reflectance between 750 and 1800 nm and about 2000 and 2400 nm. A narrow band of lower reflectance is observed at about 1900 nm. In Fig.13, spectral response of the exhaust pipe is increasing with the wavelengths. Indeed, the spectral response is low, 10% for the wavelength 400 nm while the spectral response for the same angle is 70% for the wavelength 2400 nm. Its shiny appearance confirms the reflectance peak at 10 ˚.

  Figures 8 to[START_REF] Marchetti | Retrieving visibility distance in fog combining infrared thermography, Principal Components Analysis and Partial Least-Square regression[END_REF] show from the left to the right, the object, its spectral response and the angular distribution of the reflectance. The analysis of each object reflectance gives the following results.In Fig.8, the grey part of the cone is highly reflective for an angle of 10 ˚. The reflectance does not change with the wavelength. Fig.9shows that the reflector reflectance is decreasing with wavelength. The most reflected radiation is that of the visible. Beyond 1550 nm, it reflects only 50% of the radiation, while its spectral response is negligible at 2400 nm. It is therefore necessary to pay attention to this type of equipment, especially with Lidars using a 1550 nm laser. In Fig.10, the asphalt is not very reflective and it has a "quasi-Lambertian" behavior. Its spectral response does not exceed 10% in the infrared, with a maximum of 7.5% in the visible range. For the dark blue polar sweater, infrared wavelengths between 800 nm and 1550 nm are the more reflecting (Fig.11). The tire in Fig.12has a higher reflectance between 750 and 1800 nm and about 2000 and 2400 nm. A narrow band of lower reflectance is observed at about 1900 nm. In Fig.13, spectral response of the exhaust pipe is increasing with the wavelengths. Indeed, the spectral response is low, 10% for the wavelength 400 nm while the spectral response for the same angle is 70% for the wavelength 2400 nm. Its shiny appearance confirms the reflectance peak at 10 ˚.

ACKNOWLEDGMENT

The authors would like to thank Jean-Luc Bicard and Anthony Bayle for their active contribution to experiments and data processing.

*The research leading to these results has received funding from the European Union under the H2020 ECSEL Program as part of the DENSE project, contract number 692449. The French government research program "Investissements d'Avenir" also contributed through the IMobS3 Laboratory of Excellence (ANR-10-LABX-16-01).


Marcin P Lodzień

Maciej Kościelski

Alice Sinatra

Producing and storing spin-squeezed states and Greenberger-Horne-Zeilinger states in a one-dimensional optical lattice

published or not. The documents may come

Producing and storing spin-squeezed states and Greenberger-Horne-Zeilinger states in a one-dimensional optical lattice Marcin P Lodzień, Maciej Kościelski, Emilia Witkowska, Alice Sinatra

I. INTRODUCTION

Ultra-cold atoms are widely recognized as a platform for fundamental tests of quantum mechanics, quantum information and quantum networks [1][2][3][4][5], as well as sensors in well established applications like atomic clocks [6][7][8]. In some implementations, non-trivial correlations among the atoms in the form of squeezing and entanglement are directly profitable. In addition, recent developments in experimental techniques allow to produce and control ultra-cold systems in the few-body regime with a good precision on the atom number [9][10][11][12]. These experimental developments in ultra-cold few body systems boost theoretical studies [13] indicating new possibilities for thermometry [14] and quantum simulations [15,16].

In a previous work [17] we proposed a method for the dynamical generation of spin-squeezing and its storage in an optical lattice with unit filling, see Fig. 1. The scheme makes use of interactions in a Bose-Einstein condensate with two internal states. The atoms are initially in the superfluid phase delocalized over the entire lattice and in a single internal state. As soon as the atoms are put in a superposition of two internal states by an electromagnetic pulse, atomic interactions dynamically generate non-trivial correlations and spin-squeezing [18,19]. The lattice depth is then adiabatically increased. As the system approaches the Mott-insulator phase, the spin dynamics slows down to stop completely at the transition, while the spin-squeezing at its best and the quantum correlations survive and are stored deeply in the Mott insulating phase. The work [17], that was for a large number of atoms N ∼ 10 5 and a three-dimensional lattice, was relying on two main approximations: (i) the hypothe-FIG. 1. (a) Initially, ultra-cold atoms in an internal state a are prepared in the superfluid phase in a shallow optical lattice with unit filling. (b) At time t = 0, a π/2-pulse puts each atom in a coherent superposition of two internal states a and b. Immediately after the pulse, while the system in the superfluid phase, binary atomic interactions between atoms start the generation of quantum correlations in the system. (c) The lattice height is gradually increased, in such a way that the system undergoes the Mott-insulator transition with one atom per site at the "best time", to store in the lattice either the best squeezing or a GHZ state.

sis of adiabaticity while ramping the optical lattice in the two-component interacting condensate, and (ii) the Gutzwiller approximation to calculate the ground state energy of the system for different lattice heights.

In the present paper, we deal with a few-bosons in a one-dimensional optical lattice. Within a full represen-tation of the system Hamiltonian [20], we study its exact state and dynamics. This allows us to confirm the validity of the scheme beyond the approximations we used in [17]. The advantage of the few-body system over the macroscopic number of atoms is furthermore the possibility to access and fully characterize the quantum state stored in Mott-insulator phase, including the excitation generated by the system manipulations (the initial pulse and the dynamic ramping of the lattice) and by thermal fluctuations. From the experimental point of view, the on-site and inter-site correlation functions that we calculate would be accessible in an experiment with a quantum gas microscope [21,22].

After introducing in Section II the system Hamiltonian and parameters, the following Section III is devoted to an analytical one-axis twisting model, issued from the zero-momentum part of the two-component (two internal states) Bose-Hubbard Hamiltonian, that catches the main features of the dynamics. By testing the model against exact numerical simulations, we prove its validity for a large enough condensed fraction. One useful result, predicted by the simple model and confirmed by the simulations, is that the squeezing dynamics is faster in a shallow lattice with respect to the homogeneous system. This acceleration is directly related the confinement of the atoms causing an increase of the interaction energy and a consequent reduction of the squeezing time scales.

In Section IV we study the structure of the states stored in the Mott phase that turns out to be non trivial. With respect to what we shall call an "ideal Mottinsulator phase-state", that has exactly one atom per site, each atom being in the same coherent superposition of two internal states, the stored entangled state shows additional phase factors that play a crucial role in determining its properties. We discuss the possibility of storing either a gaussian spin-squeezed state or a maximally entangled Greenberger-Horne-Zeilinger state [23,24], superposition of two ideal Mott-insulator phase-states with different phases, provided that the ramp time is properly chosen. The challenge of generating experimentally a Greenberger-Horne-Zeilinger state has been already taken-up in several platforms, see e.g. [25][26][27][28][29][30][31]. The strength of the method we propose here relies in its relative simplicity, and in the possibility to scale it up to a larger system size, by using a larger lattice and by going from one to three spatial dimensions.

Finally, in Section V, we discuss the role of a finite temperature resulting in particle-hole excitations in the Mott phase that can lead to double occupation of certain lattice sites. We show that it is nevertheless possible to maintain squeezing in the Mott phase at low temperature, when the probability for particle-hole excitation is weak. Concluding remarks and a summary are given in Section VI.

II. CONSIDERED SYSTEM AND MODEL HAMILTONIAN

In our study, we use the Bose-Hubbard model to describe a system composed of a few ultra-cold bosonic atoms in an optical lattice potential. In the following we discuss the model, parameters that we used and the protocol for dynamical generation of correlations among the atoms.

A. Bose-Hubbard Hamiltonian

We consider an ultra-cold gas composed of a few bosonic atoms with repulsive interactions in the two internal states a and b, loaded in a three-dimensional optical lattice or optical tweezers potential. We assume that tunneling is possible only along the x direction while transversely the atoms are in a localized wave-function that is approximated by a Gaussian with a characteristic length L ⊥ . The optical lattice potential in the xdirection is V (x) = V 0 sin 2 (kx), where V 0 is the lattice height, k = 2π/λ is the wave number, and d = λ/2 is the spatial period of the lattice. In the absence of the lattice, the system is homogeneous, the atoms being confined in a flat-bottom potential [32]. After reduction of the perpendicular directions, the effective one-dimensional system Hamiltonian reads

Ĥ = σ=a,b dx Ψ † σ (x) ĥ Ψσ (x) + g σσ 2 Ψ2 † σ (x) Ψ2 σ (x) + g ab dx Ψ † a (x) Ψ † b (x) Ψb (x) Ψa (x) (1)
with the single-particle Hamiltonian ĥ = -

2 2m d 2 dx 2 + V (x). (2)
The interaction coefficients g σσ ′ = g 3D σσ ′ /(2πL 2 ⊥ ) are determined by the transverse confinement length L ⊥ and by the coupling constants in three dimensions g 3D σσ ′ = 4π 2 a σσ ′ /m, where a σσ ′ is the s-wave scattering length for one atom in the state σ and one in the state σ ′ , m is the atomic mass and is the Planck constant. We assume repulsive interactions between atoms in the two states, with a ↔ b symmetry leading to g aa = g bb , and an adjustable inter-species coupling g ab , restricting to the phase-mixing regime : g ab < g aa [33].

The system (1) in the lowest energy band is conveniently considered in the basis of Wannier functions w(xx i , t) localized around the lattice sites, where x i denotes position of the i-th site [2]. In the tight-binding limit when the lattice height is larger than the recoil energy E R = 2 k 2 /(2m) and the Wannier functions are well localized around each lattice site, the tunneling and interactions terms fall off rapidly with the distance |x i -x j |, leading to the Bose-Hubbard model

ĤBH = -J i,j=i±1 â † i âj + b † i bj + U aa 2 i na i (n a i -1) + U bb 2 i nb i (n b i -1) + U ab i na i nb i , (3)
where J and U σσ ′ are the tunneling and interaction parameters. âi ( bi ) is the annihilation operator of an atom in internal state a (b) in the i-th site of the lattice, and na i = â † i âi , nb i = b † i bi are the corresponding number operators. The Bose-Hubbard Hamiltonian commutes with the total number of particles in each component ĤBH , Nσ = 0, where Nσ = i nσ,i with σ = a, b, but it does not commute with the occupation numbers na,i , nb,i of the i-th site, due to the presence of the hoping terms.

We address here the case in which the total filling is commensurate with the lattice. In the repulsive interaction regime we consider, a transition from a superfluid to a Mott-insulator phase is expected in the two-component system, as in the usual case of a single species [21]; the magnetic order associated with the pseudospin degrees of freedom (boson components) being a further feature in the ground state [34]. We briefly present the construction of the Bose-Hubbard Hamiltonian (3) from the more general one (1) in Appendix A. The basic properties of the tunneling and interaction parameters are summarized in Appendix C, where we also calculate for different lattice heights the couplings appearing in (1) that are neglected in (3).

B. Procedure and parameters used in the simulations

Throughout the paper, we will consider two situations: (i) the case of a static lattice with a time independent height V 0 , and (ii) the dynamical case where the lattice is raised with a linear ramp in a time τ

V 0 (t) = V ini + (V fin -V ini ) t τ , (4)
from an initial value V ini where the system is in the superfluid regime to a final value V fin in the Mott-insulator regime. lattice sites for atoms in the internal state a (b). According to the commutation relations ĤBH , Nσ = 0, the Hilbert space can be written as a direct sum of subspaces with a fixed number of atoms in each internal state . A numerical method for the Fock state basis generation is described in [20]. Operators Ô are represented by matrices n ′ | Ô|n , and pure states by vectors whose components are the coefficients of the decomposition in the basis.

κ = κ N,0 ⊕ κ N -1,1 ⊕ κ N -2,2 • • • κ 0,N ,
We consider a spin coherent state as the initial state for the squeezing evolution. In order to prepare it numerically, we assume that initially all the atoms are in the internal state a and, in the case of zero temperature, in the ground state of H BH for the initial values of J and U :

Ψ 0 (0 -) = {n a i } i n a i =N α {n a i } |{n a i }, {0} , (5)
where the sum runs over all the possible configurations, and the coefficients α {n a i } are calculated numerically. A π/2 pulse, equivalent to a rotation of the state around the y-axis through angle π/2, is applied to create the spin coherent state [35] |Ψ (0 + ) = e -i Ŝyπ/2 |Ψ 0 (0 -) .

When the temperature is nonzero, the initial state after rotation is given by the density matrix:

ρT (0 + ) = e -i Ŝyπ/2 ρT (0 -)e i Ŝyπ/2 , (7)
with

ρT (0 -) = 1 Z n e -En/kB T |Ψ n (0 -) Ψ n (0 -)| (8)
where E n are eigenenergies of the initial Hamiltonian with i n a i = N and |Ψ n (0 -) are the corresponding eigenvectors, k B is the Boltzman constant, T the temperature and Z the normalization. The initial state is then evolved according to the Schrödinger equation i ∂ t |Ψ (t) = ĤBH |Ψ (t) for the zero temperature case, or the von Neumann equation i ∂ t ρT (t) = [ ĤBH , ρT (t)] for the nonzero temperature case. Both equations were solved numerically using the Runge-Kutta method.

In Fig. 2 we show the values of the tunneling parameter J(t) and the interaction parameter U aa (t) through the considered ramp (4) optimized to produce and store a spin-squeezed state with N = 6 atoms. J(t) and U aa (t) were calculated numerically according to the definitions (A3a-A3b) using the Wannier functions, as detailed in Appendix A. The numerical values of parameters in all the simulations, are typical for hyperfine states in alkali atoms 1 . We took a a = 100.4a B , a b = a a , and a ab = 95a a , where a B = 5.29 × 10 -11 m is Bohr radius, the lattice constant d = λ/2 = 431nm, the ratio (a a + a b -2a ab )/d = 1.25 × 10 -3 , a a /d = 0.012. The recoil energy value is E R = 2.04 × 10 -30 J and the characteristic time scale t 0 = /E R = 51.5µs. We also assume that L ⊥ = d/ √ 2π.

C. Spin squeezing parameter

We quantify the level of spin-squeezing by the parameter [37,38]

ξ 2 = N ∆ 2 Ŝ⊥ min S 2 , ( 9
)
where S is the length of the mean collective spin and ∆ 2 Ŝ⊥ min is the minimal variance of the spin orthogonally to the mean spin direction. In the multimode case of atoms in an optical lattice, collective spin operators are defined as a sum of local spin operators

Ŝx = M i=1 1 2 â † i bi + b † i âi , (10) Ŝy
= M i=1 1 2i â † i bi -b † i âi , (11) Ŝz
= M i=1 1 2 â † i âi -b † i bi (12)
As for individual spins, collective spin components obey cyclic commutation relations Ŝx , Ŝy = i Ŝz . Due to the 1 In an experiment, one could use for example 87 Rb atoms in states |F = 1, m F = 1 and |F = 2, m F = -1 where the interspecies scattering length a ab is tuned by means of a Feshbach resonance [19], or the states

|F = 1, m F = -1 and |F = 2, m F = -2
where the interspecies interaction can be tuned by slightly shifting the optical lattices for the two components [36].

fact that the Bose-Hubbard Hamiltonian commutes with Ŝz and the initial state is an eigenstate of Ŝx one can write explicitly the expression for minimal fluctuations

∆ Ŝ2 ⊥,min = 1 2 2 ∆ Ŝ2 z + A -A 2 + B 2 , (13)
A = ∆ Ŝ2 y -∆ Ŝ2 z , (14)
B = 2Re Ŝy Ŝz -Ŝy Ŝz . (15)

III. CONSTRUCTION AND VALIDITY OF TWO-MODE MODEL

In this section we introduce a two-mode model that can be used to obtain analytical predictions [17], and investigate the dependence of the best squeezing time on the lattice potential height. To construct the two-mode model in the most general way, we introduce the ground state energy E 0 ( Na , Nb ) of the system having N a atoms in the internal state a and N b atoms in b, and make its Taylor expansion up to the second order around the average values Na and Nb of the atom numbers, in the initial state at t = 0 + . This leads to

E 0 ( Na , Nb ) = E 0 + ∂ a E 0 ( Na -Na ) + ∂ b E 0 ( Nb -Nb )+ + 1 2 ∂ 2 aa E 0 ( Na -Na ) 2 + 1 2 ∂ 2 bb E 0 ( Nb -Nb ) 2 + + 1 2 (∂ 2 ab E 0 + ∂ 2 ba E 0 ) ( Na -Na )( Nb -Nb ) + • • • , (16)
where

E 0 ≡ E 0 ( Na , Nb ), ∂ σ E 0 ≡ ∂E0 ∂Nσ | Na, Nb and ∂ 2 σσ ′ E 0 ≡ ∂ 2 E0 ∂Nσ∂N σ ′ | Na, Nb , with σ, σ ′ = a, b.
Notice, that the first derivative of the ground state energy is the zero temperature chemical potential µ σ = ∂ σ E 0 of the σ component. By introducing the collective spin components, the ground state energy expanded to the second order ( 16) becomes [39]

E 0 ( Na , Nb ) = f N + v N Ŝz + χ Ŝ2 z (17)
In the above equation, f N is a function of the total number of particles, the linear term in Ŝz describes spin precession around the z axis with the velocity v N , and finally

χ = 1 2 (∂ a µ a + ∂ b µ b -∂ a µ b -∂ b µ a ) . (18)
In the symmetric case that we consider, with U aa = U bb , one has v N = 0 (no precessions). Note that the approximation of the energy by its second order expansion is more and more accurate as N is large, as the relative width of the distributions of N a and N b decrease as 1/ √ N . The approximated energy (17) is the basis of the twomode model we consider here. It has the form of the oneaxis twisting (OAT) model introduced by Kitagawa and Ueda [40] and it catches the main features of the squeezing dynamics. In particular we expect that the maximum level of squeezing generated in our scheme is the same as the one generated by equivalent the OAT model (17), in both the static and the dynamic scheme that we consider. The time scale, or the best squeezing time, depends on the system Hamiltonian trough the parameter χ. Here below, we give an approximate expression of χ in the static case, assuming that all the atoms are in the zero momentum mode that is the condensate mode in the noninteracting case. To this aim, we write the Bose-Hubbard Hamiltonian (3) in the momentum representation, âi = 1 √ M q e -iqxi âq and bi = 1

√ M q e -iqxi
bq , and we keep only the zero momentum terms [41]. We obtain :

ĤBH, 0 ( N0,a , N0,b ) = ǫ 0,a N0,a + ǫ 0,b N0,b + U aa 2M N 2 0,a + U bb 2M N 2 0,b + U ab M N0,a N0,b , (19)
where

N0,a = â † q=0 âq=0 , N0,b = b † q=0 bq=0 , ǫ 0,a = -2J - U aa /(2M ) and ǫ 0,b = -2J -U bb /(2M ). The chemical potentials are µ a | N0,a, N0,b = ǫ 0,a + Uaa M N0,a + U ab M N0,b , µ b | N0,a, N0,b = ǫ 0,b + U bb M N0,a + U ab M N0
,a where all expressions are derived assuming that all atoms are in the zero momentum mode. One then obtains

χ q=0 = 1 2M (U aa + U bb -2U ab ), (20)
Using (20), the parameter χ q=0 can be calculated numerically using the exact form of the Wannier function to determine the interaction parameters U σσ ′ . Moreover, analytical expressions for χ q=0 can be obtained in the limiting cases of a very shallow or a very deep lattice. The first limit is that of a homogeneous system when V 0 → 0. The values of the interaction parameters are then U homo σ = g σ N/L x , where L x = M d is the system size along x, leading to

χ homo q=0 = 4 M π ∆a d E R (21)
where ∆a = a a + a b -2a ab . The second one is that of localized Wannier functions that can be approximated by Gaussians, see Appendix C. In this case one obtains

χ Gauss q=0 = 32 π ∆a d V 0 E R 1/4 E R (22)
The zero momentum part of the Bose-Hubbard Hamiltonian (19) can be viewed as bimodal: the atoms in q = 0 of the a state and those in q = 0 of the b state. The evolution of such bimodal subsystem can be solved analytically in the Fock state basis: The black solid lines show the exact numerical results while the red dashed lines the prediction of 2MM with χq=0 from (20). The agreement between the 2MM and exact results degrades for larger V0. taking the spin coherent state as the initial state that determines the values of c N0,a = 2 -N/2 N N0,a . In (23) we omitted the constant global phase factor. Consistently, for the two-mode model, the spin squeezing parameter ( 9) is considered for the zero momentum mode, replacing the collective spin operators Ŝx Ŝy and Ŝz by their zero momentum counterparts Ŝ0,x = (â † q=0 bq=0 + b † q=0 âq=0 )/2, Ŝ0,y = (â † q=0 bq=0 -b † q=0 âq=0 )/2i and Ŝ0,z = ( N0,a -N0,b )/2. The evolution of such spin squeezing parameter for zero-momentum-mode is solvable analytically [40,42]. We will refer to it as the two-mode model (2MM) while comparing it to the exact numerical results. We expect agreement between the exact calculation and the 2MM as long as N0,a + N0,b ≈ N in the initial state before π/2 pulse and as long as the system stays in the superfluid regime.

|ψ 0 (t) = N N0,a=0 c N0,a e -i χ q=0 t 4 (N -2N0,a) 2 |N 0,a , N -N 0,a , (23)
In Fig. 3 we show an example of the evolution of the spin squeezing parameter (9) for two different static values of the lattice height V 0 /E R = 0.4 and V 0 /E R = 6, both for N = 6. The maximum level of squeezing achieved during the evolution is the same for the two values of V 0 /E R while the time scales are different, as expected. The agreement between the 2MM and the exact results degrades when V 0 increases, because the approximation of having all the atoms in the zero momentum mode, that we used to calculate χ q=0 , becomes unjustified. To be more specific on this point, in Fig. 4 we collect the best squeezing times calculated exactly (black solid line) and compare them to the 2MM results (orange dot-dashed line) with χ q=0 calculated from (20). The condensed fraction defined as the average occupation of the q = 0 momentum mode:

f σ c = 1 N σ † q=0 σq=0 = 1 N M i,j σ † i σj , (24)
where σ = a, b is shown on the same figure. The following observations can be made: (i) the 2MM works well as long as the condensed fraction is close to one, and (ii) the best squeezing time is smaller than for an homogeneous system with the same number of atoms in a specific range of V 0 . The squeezing acceleration by the optical lattice potential in the range 0 4. The upper bound for V 0 expresses the fact that the lattice height should be low enough for the system to remain in the superfluid phase with a condensed fraction close to one. To obtain an estimate of that for small atom numbers, we can use the explicit analytical expression of the condensed fraction for

.4E R < V 0 < 7E R ,
N = 2, see Appendix D. It is f a c = 64J 2 +16JΩaa++Ω 2 aa+ 2(64J 2 +Ω 2 aa+ ) with Ω aa+ = U aa + 64J 2 + U 2 aa . It is obvious that f a c → 1 when U aa /J → 0.
In order to estimate the validity range of the 2MM we can fix the condensed fraction to a value close to one. By taking for example f a c ≈ 0.99, the threshold point above which the 2MM fails is then (U aa /J) th ≈ 1.6 corresponding to (V 0 /E R ) th ≈ 7, which gives a fair estimate of what we see in the exact numerical calculations in Fig. 4.

IV. SPIN-SQUEEZED AND GHZ STATES STORED IN MOTT PHASE

In order to store the best squeezing in the Mottinsulator phase, we linearly ramp the lattice potential depth according to (4). We choose the initial value of the lattice depth V 0 = 3E R as it corresponds to a condensed fraction close to unity, and it generates the squeezing faster than for a homogeneous system. An example of spin squeezing parameter and condensed fraction evolution 2 is shown in Fig. 5. The squeezing parameter stops evolving around its minimum because the ramp is adjusted so that the system undergoes the Mott transition at the best squeezing time. In the Mott-insulator phase with one atom per site, atoms are isolated around lattice sites and stop to interact. As we will show it in this section, not only a squeezed state can be stored in the Mott phase. By varying the ramp time τ one can store other states dynamically produced by the interactions, including macroscopic superpositions of phase states.

To gain understanding about the structure of state stored in Mott phase, lets us first consider the simplest case of N = 2 and then generalize the arguments for an arbitrary number of atoms. 2 The values of the tunneling and interaction terms were calculated exactly according to (A3a) and to (A3b). They are shown in In general, the state of the system can be considered in the Fock state basis as we did numerically. For N = 2 it reads

|Ψ(t) = 10 w=1 c w (t)|w , (25)
The coefficients of the decomposition in the Fock basis satisfy c 6 = c 5 , c 7 = c 4 , c 8 = c 3 , c 9 = c 2 , c 10 = c 1 , c 3 = c 1 and c 8 = c 10 , as a consequence of the symmetry a ↔ b with respect to the exchange of the internal states and of the symmetry with respect to exchange of the lattice sites. The evolution of the coefficients can be found analytically when U aa and J are time-independent as discussed in Appendix D. In the time-dependent case, the evolution of the c w can be treated numerically and the result is shown in Fig. 6. It turns out that the state stored in the Mott-phase has an interesting structure. The coefficients c w are nonzero only for Fock states that have one atom per lattice site. All the nonzero coefficients have identical absolute values and the only difference among them is in the phase factor. Referring to the states (26), we find numerically that the state stored in the Mott phase at the end of the ramp for N = 2 has the form

|Ψ MI,N =2 = 1 2 |2 + e iφ |5 + e iφ |6 + |9 (27)
up to a constant global phase factor. The only nonzero coefficients in the long time limit are c 2 , c 5 , c 6 and c 9 .

The very same Fock states appear in what we shall call the ideal φ = 0 Mott-insulator phase-state, 28) obtained from a single component Mott state, by putting each atom in the coherent-superposition (|a + |b )/ √ 2. However, unlike in the φ = 0 Mott-insulator phase-state, a nonzero relative phase appears between individual Fock states. The value of the relative phase φ depends on the state that is stored in the Mott phase. The above introduced structure (27) in term of Fock states, is illustrated in Fig. 7 for different values of the ramp time, corresponding to different quantum states stored in the Mott-insulator phase.

|φ = 0 MI = a † 1 + b † 1 √ 2 a † 2 + b † 2 √ 2 |0 = |2 + |9 + |5 + |6 2 (
The general form (27) of the stored state for N = 2 can be expressed in terms of two-site phase-states,

|φ 1 , φ 2 = a † 1 + e iφ1 b † 1 √ 2 a † 2 + e iφ2 b † 2 √ 2 |0 . (29)
Indeed one has

|Ψ MI,N =2 = e iφ/2 cos φ 2 |0, 0 + isin φ 2 |π, π . (30)
It is true for any value of the phase φ. In this sense, the state stored in the Mott phase is a superposition of two different Mott-insulator phase-states with different phases, and it becomes an even superposition when φ = π/2.

B. Case N = 6

The natural question arises how the discussed above structure changes when the number of atoms increases. One might expect that the state stored in the Mott phase is a superposition of 2 N Fock states with single occupation per lattice site and N/2 different phase factors. We have checked numerically that this is the case for N = 4 and N = 6. In the latter case one has

|Ψ MI,N =6 = 1 2 3 S 6,0 + S 0,6 + e iφ51 (S 5,1 + S 1,5 ) +e iφ42 (S 4,2 + S 2,4 ) + e iφ33 S 3,3 ] , (31)
were we introduced notation |{n ia }, {n ib } denotes the multi-site Fock state. For example S 6,0 = |111111, 000000 . The structure (31) of the state for N = 6, is illustrated in Fig. 8 for different values of the ramp times, corresponding to different quantum states stored in the Mott-insulator phase.

S Na,N b = {n ia } ′ i n ia =Na {n ib } ′ i n ib =N b |{n ia }, {n ib } , (32)
As we did it previously, we can express the stored state in terms of Mott-insulator multi-site phase-states. In general, in the case of six atoms and six lattice sites, the state (31) can be expressed in terms of six-site phasestates

|φ 1 , φ 2 , φ 3 , φ 4 , φ 5 , φ 6 = 6 j=1 â † j + e iφj b † j √ 2 |0 = 1 2 3 6 Na=0 SNa,N-Na , (33)
with

SNa,N b = {n ia } ′ i n ia =Na {n ib } ′ i n ib =N b e i i φin bi |{n ia }, {n ib } .
(34) When considering a Mott-insulator phase-state |φ MI (33) where the phase is the same for all lattice sites,

|φ MI ≡ |φ, φ, φ, φ, φ, φ (35
+ | π 3 MI + | 2π 3 MI + | 4π 3 MI + | 5π 3 MI , 1 2 (S 51 + S 15 ) = |0 MI -|π MI -i | π 2 MI + | 3π 2 MI , 3 4 (S 42 + S 24 ) = 2|0 MI + 2|π MI -| π 3 MI + | 2π 3 MI + | 4π 3 MI + | 5π 3 MI , 1 2 S 33 = |0 MI -|π MI + i | π 2 MI + | 3π 2 MI . ( 36
)
In the special case when e iφ42 = 1 and φ 51 = φ 33 = π/2, one obtains the N = 6 Greenberger-Horne-Zeilinger state

|Ψ GHZ,N =6 = e iπ/4 √ 2 [|0 MI + i|π MI ] . (37)
As shown in Fig. 8, this occurs when the ramp time τ equals half of period of the one-axis twisting dynamics, corresponding to freezing in the Mott state the correlations of a Schrödinger cat, macroscopic superposition of two phase-states. Note that ( 37) is equivalent to the most common representation of the Greenberger-Horne-Zeilinger state in the number operator basis. To see this, one should change the basis by rotation of the state (37) by π/2 around the y-axis.

C. On-site and inter-site correlation functions in the Mott-squeezed state

The stored Mott-squeezed state is neither a φ = 0 Mott-insulator phase-state nor a two-mode squeezed state. It is a multimode entangled state possessing some of the properties the two. To have more insight into the character of the Mott-squeezed state, we analyze here the two correlation functions: â † i bj and { Ŝy,i , Ŝz,j } -2 Ŝy,i Ŝz,j , where Ŝy,i = (

â † i bi -b † i âi )/2i and Ŝz,i = (â † i âi -b † i bi )/2.
The first one quantifies the correlation between the internal states a and b of the atoms, that is maximal for i = j in a φ = 0 Mottinsulator phase-state. The second one characterizes spinspin correlations between different atoms, that is a property of spin-squeezed states as can be seen by the definitions ( 9) and ( 15). The higher the value of the spinspin correlation function the higher the level of squeezing. Note that Ŝy,i = Ŝz,i = 0 both for the initial state and during evolution.

In Fig. 9 we show both correlation functions for different distances ij between the sites, both in the case of a linear ramp of the lattice (solid lines) and in the case of a static lattice (dashed lines). In the case of a ramp, in the long time limit when J → 0, the | â † i bj | 2 correlation function is nonzero only on-site; the off-site correlation function â † i bj tends to zero independently of the value ij as long as i = j. On the other hand, the spin-spin correlation function is nonzero only between different sites i = j.

The origin of such a behavior comes from the structure of the Mott-squeezed state, and can be understood in detail in the case N = 2. The correlation functions of interest are in this case

â † 1 b1 = 2 √ 2|c 1 c 4 |cos(θ 1 -θ 4 ) + 2|c 2 c 5 |cos(θ 2 -θ 5 ), (38) â
† 1 b2 = 2 √ 2|c 1 c 5 |cos(θ 1 -θ 5 ) + 2|c 2 c 4 |cos(θ 2 -θ 4 ), (39)
and

{ Ŝy,1 , Ŝz,1 } = 2 √ 2|c 1 c 4 |sin(θ 4 -θ 1 ), (40)
{ Ŝy,1 , Ŝz,2 } = 2|c 2 c 5 |sin(θ 5 -θ 2 ), (41)
where we used (25), the a ↔ b symmetry and the notation c w = |c w |e iθw . Since with one atom per site, only the coefficients c 2 , c 5 , c 6 , c 9 are nonzero, the analytical expressions ( 38)- (39) confirm that â † i bj is nonzero when 0 1000 2000 3000 4000 5000 tE R / / h 0 0,1 i = j and zero for i = j. One can check that the product of nonzero coefficients (those of Fock states with one atom per lattice site) always contribute to the on-site function â † i bi , independently of the total atom number as long as unit filling is considered. Similarly one can see from (40) and ( 41) that the spin-spin correlation function has the opposite behavior, displaying a nonzero product of coefficients when i = j only3 .

V. EFFECT OF NONZERO TEMPERATURE ON THE MOTT-SQUEEZED STATE

Two types of excitations exist in the Mott phase of our two component system: gapped excitations corresponding to double occupations of a site, for which the energy scale is the on-site interaction energy U σσ or U σσ ′ , and "soft" excitations within the single occupation manifold that can be described by an effective spin model. In our protocol it is important to maintain the system in the single occupation manifold in order to keep the stored squeezing and entanglement constant in time. Within this manifold, it is not important to remain in the ground state because the spin-dependent interaction energy is strongly reduced, scaling as J 2 /U σσ ′ [34,[43][44][45].

We therefore concentrate on excitations out of the low-energy manifold toward states with double occupations. We show that thermal fluctuations cause the squeezing parameter to oscillates in the Mott phase. Its value nevertheless remains small in the low temperature limit. Residual interaction in the Mott phase due to double occupation of a site, directly influences our scheme because the squeezing dynamics is then imperfectly stopped. To give a rough estimate, we expect that thermal effects become visible when the temperature is comparable with the excitation energy for double site occupation k B T ≥ U ab ≃ U bb = U aa , where k B is the Boltzmann constant.

According to Fig. 2 and Fig. 5, we see that U aa varies between 0.065E R and 0.08E R in the critical region from tE R / = 1000 and tE R / = 2000, after which the thermally induced double occupations are frozen in the deep Mott phase. We then expect notable thermal effects for temperatures of this order. The finite temperature simulation results are shown in Fig. 10(a) for N = 4. We observe that nonzero temperature not only increases the value of the best squeezing stored in the Mott phase as it was already understood for bimodal condensate [46]; in addition, an evolution of the squeezing parameter stored in the Mott phase takes place in the form of regular oscillations. We expect the amplitude of oscillations of ξ 2 (t) in the Mott phase to be proportional to the double occupation probability, quantified by the second order correlation function g

(2) σ,σ ′ = σ † i σ′ † i σi σ′ i / σ † i σi σ′ † i σ′ i with
σ, σ ′ = a, b, illustrated in Fig. 10(b). Indeed, while the temperature changes from 0.06E R to 0.08E R , both the stationary value of the second order correlation function and the oscillation amplitude of the squeezing parameters increase by approximately the same factor. On the other hand, the period of the oscillations in ξ 2 (t) in the Mott phase is determined by the periodicity of the dynamics for two atoms on the same site, similarly to [47], corresponding here to ϕ = (U aa -U ab )t/ = 2π. Using the value of U aa at the end of ramp, for V 0 /E R = 40 one finds a value of the period of approximately 1369 /E R 4 in agreement to what observed in Fig. 10.

Our results demonstrate the possibility to maintain spin-squeezing in the Mott phase in our small system even at finite temperature. 4 At the end of the ramp, U bb = Uaa/E R ≈ 0.09, and U ab = 0.95Uaa. One then finds ∆t = 2π /(Uaa -U ab ) = 1369.22E R / . As long as the initial condensed fraction is large, and the initial temperature is week with respect to the value of U ab and Uaa (as discussed in the text), the effect of the temperature is very weak.

VI. SUMMARY

We investigated the dynamical generation and storage of spin-squeezed and Greenberger-Horne-Zeilinger states in the Mott phase of system of a cold atoms in an optical lattice. This is done by raising an optical the lattice potential in an interacting bimodal Bose-Einstein condensate that is initially in the superfluid phase in the lattice.

We introduced and explained a two-mode one-axis twisting model which allows analytical estimates, and used it to show that the presence of the optical lattice can accelerate the entanglement dynamics as compared to an homogeneous system. The careful analysis of the stored entangled state in the Mott phase revealed its structure: it is composed by the same Fock states that are present in an ideal φ = 0 Mott-insulator phase-state, but with additional relative phase factors depending on the state that is stored. We show that by a proper choice of the lattice ramp time, that would cor-respond to formation of a Schrödiger cat in a two-mode one-axis twisting model, it is possible to store in the lattice a Greenberger-Horne-Zeilinger state. The structure of the Mott-squeezed states was used to explain the behavior of several correlation functions, on-site and between different sites, that we have calculated numerically. The first-order internal state correlation functions â † i bj are nonzero only when i = j demonstrating one-body coherence within each atom. The spin-spin correlation functions { Ŝy,i , Ŝz,j } , that characterize spin-squeezing are nonzero only when i = j as a consequence of unit filling. Finally, we discussed the effect on the squeezing of nonzero temperature that allows particle-hole excitations. We show that it is still possible to store a squeezed state in the Mott phase, although oscillations of the squeezing parameter appear at finite temperature due to the nonzero probability of double occupation of a single site.

From the experimental point of view, the system with a few atoms is very interesting because of the large control one has on the system, both for the preparation and for the measurement. A step-by-step scaling up of the system size, by increasing the number of lattice sites or the spatial dimension, would allow to explore the boundary between the microscopic and the macroscopic world.

+ V 0 2E R a j - V 0 4E R (a j+1 + a j-1 ) = E E R a j .
(B5) This is a standard eigenvalue problem and a numerical solution of the polynomial equation det(M -E1) = 0 with

M j,j ′ E R =    2j + q k 2 + V0 2ER , if j = j ′ ; -V0 4ER , if |j -j ′ | = 1; 0, otherwise, (B6)
provides (2j max + 1) eigenvalues of the matrix M . These eigenvalues are labeled by the index l = 1, 2, ..., 2j max + 1 (called band index) and denoted as E l,q . The components of an eigenvector corresponding to given eigenvalue E l,q are the coefficients a j of (B3). Having them, one needs to evaluate the sum in (B3) to extract the Bloch function (B2).

Appendix C: Properties of tunneling and interaction terms

In the limit V 0 /E R ≫ 1 the Wannier functions can be approximated by Gaussian functions whose width is set by the frequency associated to the each lattice site minimum [2]. Therefore, w(x) ≈

k 2 π 1/4 V0 ER 1/8 e - √ V0k 2 x 2 /2
. This is a fairly good approximation to obtain the interaction coefficient :

U Gauss σσ ′ (t) E R ≈ 32 π a σσ ′ d L x L y V 0 E R 1/4 . (C1)
In Fig. 11(a) we show the interacting terms calculated exactly from (A3b) and compare them to the approximated formula (C1). One can easily see that the interaction terms beyond the terms involving nearest neighbors can be neglected as compared to the latter ones, when V 0 /E R > 1 as expected [2].

In the case of the hopping parameters J(ij), the Gaussian approximation gives a relative error growing with the lattice height [48]. The hopping matrix element is rather conveniently approximated by the width of the lowest band in the one-dimensional Mathieu equation [49], yielding to

J app E R ≈ 4 √ π V 0,x E R 3/4 e -2 V 0,x E R , (C2)
for nearest neighbors, i.e. ij = 1. In Fig. 11(b) we present the tunneling terms calculated exactly and using the above formula. Indeed, as long as V 0 /E R > 1 the leading nearest neighbor term dominates. (D3)

In order to calculate the initial spin-coherent state, let us first consider the ground state when all atoms are in the a state, and then rotate it through π/2 around the y-axis of the Bloch sphere. When all atoms are in the component a then it is enough to consider the 3×3 matrix (D2). Eigenvalues of that matrix are E 1 = Ω aa-/2, E 2 = U aa , E 3 = Ω aa+ /2, where Ω aa± = U aa ± 64J 2 + U 2 aa , with the corresponding states:

|Ψ 1 = 1 N 1 |1 + Ω aa+ 4 √ 2J |2 + |3 (D4) |Ψ 2 = 1 √ 2 (-|1 + |3 (D5) |Ψ 3 = 1 N 3 |1 + Ω aa- 4 √ 2J |2 + |3 , ( D6
)
where N 2 1 = 2+Ω 2 aa+ /(32J 2 ) and N 2 3 = 2+Ω 2 aa-/(32J 2 ). The ground state before rotation is therefore 2J , 1, 0, 0, 0, 0, 0, 0, 0) T /N 1 , as E 1 has the lowest value.

GS -= |Ψ 1 =
The π/2 pulse e -i Ŝyπ/2 in the basis we consider is when the coefficients in the Bose-Hubbard Hamiltonian are time-independent, or has to be calculated numerically in the time-dependent case. The time-independent case (D11) can be solved analytically, however the general expressions are quite complex. Here, we give a result for the symmetric case when U aa = U bb , which is and c 3 (t) = c 8 (t) = c 10 (t) = c 1 (t), c 6 (t) = c 5 (t), c 7 (t) = c 4 (t), c 9 (t) = c 2 (t). In the above expressions ω σσ ′ = √ 64J 2 + U σσ ′ and Ω σσ ′ ± = U σσ ′ ± ω σσ ′ . The reason why coefficients c 6c 10 are expressed in terms of c 1c 5 is the symmetry in respect to exchange of a and b, forced by equality of interaction coefficients U aa = U bb . Otherwise, they are different and have complex analytical forms. Note, in the limit J → 0 coefficients c 2 , c 5 , c 6 , c 9 are nonzero while the remaining ones c 1 , c 3 , c 4 , c 7 , c 8 , c 10 tend to zero. Therefore, in the Mott phase only the Fock states with one atom per lattice site contributes to the nonzero values of an observable. The same holds in the non symmetric case as well.

                  
|Ψ(t) = w c w (t)

  where κ Na,N b denotes the subspace with N a atoms in the state a and N b in the state b. The dimension d Na,N b of each subspace depends on N a , N b = N -N a and on the number of lattice sites M . One has κ Na,N b = Na+M-1 Na N b +M-1 N b

FIG. 3 .

 3 FIG. 3. Squeezing parameter (9) versus time for N = M = 6. (a) V0 = 0.4ER and fc > 0.99. (b) V0 = 6ER and fc = 0.97.The black solid lines show the exact numerical results while the red dashed lines the prediction of 2MM with χq=0 from(20). The agreement between the 2MM and exact results degrades for larger V0.

FIG. 4 .

 4 FIG.4. Best squeezing time versus the static value of V0/ER extracted from the exact numerical simulations (black solid line), from the 2MM (orange dot-dashed line) and for a homogeneous system (green dashed line). The corresponding condensed fraction(24) in state a before the π/2 pulse (red solid line) is also shown, with the scale on the right axis. The vertical black dotted line marks V0/ER = 0.4 where the presence of a lattice starts to accelerate the squeezing dynamics.

FIG. 5 .

 5 FIG. 5. (a) Squeezing parameter ξ 2 and (b) condensed fraction f a c , for N = 6, are shown (i) across the lattice ramp (4) with τ ER/ = 4255 (black solid lines), and (ii) in a static situation for V0/ER = 3 (green dashed lines). The same quantities in the π/2 "Mott-insulator phase-state" (28), for which ξ 2 MI π/2 = 1 and f a c MI π/2 = 1/N , are shown for comparison (blue dot-dashed lines).

2 FIG. 6 .

 26 FIG. 6. Time evolution of the coefficients in the state (25), for N = 2 and τ ER/ = 4255. The moduli |cw(t)| are in the main figure while the inset shows the relative phases of the cw(t) that are nonzero in the long time limit. Here we introduced the notation cw = |cw|e iθw .

where|1=

  |20, 00 , |2 = |11, 00 , |3 = |02, 00 , |4 = |10, 10 , |5 = |10, 01 , |6 = |01, 10 , |7 = |01, 01 , |8 = |00, 20 , |9 = |00, 11 , |10 = |00, 02 .

FIG. 7 .

 7 FIG.7. N = 2. (a) Squeezing parameter versus ramp time τ , (b) corresponding amplitudes and (c) relative phases divided by 2π of nonzero coefficients of the state decomposition in the Fock basis, in the long time limit after stabilization in the Mott phase. The value of the squeezing parameter stored depends on the ramp time showing that an appropriate choice of τ is crucial to efficiency of the protocol. The dashed black vertical line marks the best squeezing. The only nonzero coefficients are c2, c5, c6 and c9 and the nonzero relative phases : θ5 -θ2 and θ6 -θ2 are equal, demonstrating validity of the state(27). This relative phase, marked by φ in(27), linearly depends on τ . The vertical dot-dashed orange line marks the ramp time τ corresponding to φ = π/2 for storing the N = 2 even superposition state(30).

FIG. 8 . 4 (

 84 FIG.8. N = 6. (a) Squeezing parameter versus ramp time τ , (b) corresponding amplitudes and (c) relative phases divided by 2π of nonzero coefficients of the state decomposition in the Fock basis, in the long time limit after stabilization in the Mott phase. There are 64 nonzero coefficients cw whose absolute values are close to 2 -3 . The three different relative phases between nonzero coefficients, φ33, φ42, φ51 in(31), linearly depend on τ with individual slopes. The ramp times τ for storing the best squeezed and the GHZ state are marked by the vertical black dashed and orange dot-dashed lines, respectively.

FIG. 9 .

 9 FIG. 9. (a) Internal states correlation functions | â †i bj | 2 and (b) spin-spin correlation functions { Ŝy,i, Ŝz,j} , for the linear ramp (solid lines). For comparison, the same correlation functions are also shown for the static situation corresponding to the initial value of V0/ER = 3 (dashed lines). Correlation functions between different sites i -j are marked by different colors as indicated in the legends.

FIG. 10 .

 10 FIG. 10. (a) Squeezing parameter versus time at nonzero temperature and (b) second order correlation function. Here N = 4. The values of the initial condensed fractions are f a c = 0.997 (kBT /ER = 0.02), f a c = 0.996 (kBT /ER = 0.04), f a c = 0.989 (kBT /ER = 0.06) and f a c = 0.973 (kBT /ER = 0.08).As long as the initial condensed fraction is large, and the initial temperature is week with respect to the value of U ab and Uaa (as discussed in the text), the effect of the temperature is very weak.

FIG. 11 .

 11 FIG. 11. (a) Interaction terms U aa i,j,k,l / [a σσ ′ d/(LxLy)] defined in Eq. (A3b) calculated numerically using the Wannier function for the values of i, j, k, l as indicated in legend (solid lines). The black dashed line shows the result of Gaussian approximation for U0000/ [a σσ ′ d/(LxLy)] (C1). (b) Tunneling terms in units of the recoil energy ER calculated exactly from the energy spectrum of the single-particle Hamiltonian in the case of different i -j sites distances as indicated in legend. The approximated result for the nearest neighbor (C2) is shown by the orange dot-dashed line.

Appendix D :

 : Evolution for N = 2: analytical results below we present the semi-analytical method used to describe the dynamics of N = M = 2 system. Let us consider the Bose-Hubbard Hamiltonian (3) in the Fock state basis, which for N = M = 2 is ordered and named as follows |1 = |20, 00 , |2 = |11, 00 , |3 = |02, 00 , |4 = |10, 10 , |5 = |10, 01 , |6 = |01, 10 , |7 = |01, 01 , |8 = |00, 20 , |9 = |00, 11 , |10 = |00, 02 . in the above basis has block diagonal form ĤBH = diag{ Ĥaa , Ĥab , Ĥbb } where matrix representation of Ĥaa and Ĥbb in states {|1 , |2 , |3 } and {|8 , |9 , |10 }, respectively, has form Ĥaa representation in states {|4 , |5 , |6 ,

  c -= (1, Ωaa= 4 √

  leading to the spin coherent state|Ψ(0) = e -i Ŝyπ/2 GS -= w c w (0) |w ,(D9)with c w (0) =

  is simply explained in the following paragraph. Lets us first consider the condition t best < t homo best . It implies χ homo < χ and hence λ 2 |w| 4 dx > 1. By approximating the Wannier function w by a Gaussian one obtains the lower bound for the lattice height V 0 /E R > ( 2

π ) 2 ≈ 0.4. This number is in a good agreement with what we observe in the numerical calculation, see vertical dotted line in Fig.

  Ωaa+ J . A solution of the Schrödinger equation i ∂ t |Ψ(t) = H |Ψ(t) reads |Ψ(t) = e -itH/ |Ψ(0) ,

	(2	√ 2,	a 2	, 2	√ 2, 4,	a 2	,	a 2	, 4, 2	√ 2,	a 2	, 2	√ 2) T	1 8 2 + a 2 √ (D10)
	and a = (D11)

  (t) = 2Je -itΩ ab+ /2 Ω ab+ -Ω aa+ + e itω ab (Ω aa+ -Ω ab-) + U aa Ω aa+ ) × × Ω aa+ ω ab cos(tω ab /2) + i(64J 2 + Ω aa+ U ab )sin(tω ab /2) ,

				|w ,	(D12)
	with
	c 1 (t) =	2Je -itΩaa-/2 64J 2 + U aa Ω aa+	,
	c 2 (t) =	Ω aa+ e -itΩaa-/2 2 √ 2ω aa	,
	c 4 ω 2 ab (64J 2 + Ω 2 aa+ )	,
	c 5 (t) =	e -itU ab /2 ab (64J 2 (D13) 2 2ω 2

From the mathematical point of view, it is because the two nonzero terms in the anti-commutator in(40) and (41) compensate each other for i = j. In opposite, for the j = j case they both contribute to the function as they have the same sign.

ACKNOWLEDGMENTS

We thank D. Kajtoch for very initial contribution to the project.

This work was supported by the Polish National Science Center Grants DEC-2015/18/E/ST2/00760 (MK) and under QuantERA programme, Grant No. UMO-2019/32/Z/ST2/00016 (EW). This project has received funding from the QuantERA Programme under the acronym MAQS. A.S. acknowledges support from the Horizon 2020 European Project macQsimal. M.P. was supported by the Foundation for Polish Science through the IRA Programme co-financed by EU within SG OP.

Appendix A: Bose-Hubbard model

The motion of a particle in a periodic potential formed by an optical lattice is conveniently described by the band theory, and the system Hamiltonian (1) can be considered in the basis of Bloch functions ψ l,q (x, t) [2]. These are eigenfunctions of the single-particle Hamiltonian (2) and posses required transnational properties. The Bloch functions can be constructed numerically in the plane wave basis as described in Appendix B. As a next step, the Hamiltonian is rewritten in the basis of Wannier functions w(xx i , t) localized around lattice sites, where x i denotes position of the i-th site in the lowest energy band. The Wannier functions are conveniently constructed from the Bloch states ψ l=1,q (x, t) in the following way

q∈BZ dq e -iqxi ψ l=1,q (x, t). (A1) Summation in Eq. (A1) extends over wave vectors belonging to the 1st Brillouin zone (BZ), -k < q ≤ k. The two-component Hamiltonian is then

where â † i ( b † i ) creates a particle in the single-particle Wannier state w(xx i , t) of the lowest energy band (l = 1) localized on the i-th site, in the internal state a (b). The Bose-Hubbard model considers only states in the lowest energy band, which is justified as long as the excitations energies to the higher bands are much larger than energies involved in the system dynamics. In general, if V 0 (t) is varied in time, the Wannier states, and hence the hoping and interaction parameters depend on time

q∈BZ dq E q e -i(i-j)dq (A3a)

were we introduced w(i) = w(xx i ). In Appendix C we discuss how the different tunneling (A3a) and interaction (A3b) terms depend on the lattice height V 0 .

In the tight-binding limit when the lattice height is larger than the recoil energy E R = 2 k 2 /(2m) and the Wannier functions are well localized around each lattice site, the tunneling and interactions terms fall-off rapidly with the distance |x ix j |. By keeping only the leading terms, on obtains the Bose-Hubbard model

where J = J(1), U σσ ′ = U σσ ′ 0,0,0,0 and na i = â † i âi , nb i = b † i bi . In the shallow lattice however, for V 0 ≪ 1, all the terms in (A2) have to be taken into account.

Appendix B: Numerical calculation of the Bloch functions

The eigenfunctions of the single-particle Hamiltonian (2), ĥ

are Bloch functions ψ l,q (x). They can be calculated numerically in a quite straightforward way. One starts with the representation ψ l,q (x) = e iqx u l,q (x) (B2)

where q is the quasimomentum, here an integer multiple of 2π/(M d), l denotes the energy band and x is a position. Since the functions u q (x) are periodic with the same periodicity of the lattice, they can be written as a discrete Fourier sum

while the potential V (x) = V 0 sin 2 (kx) as

(B4) Then, the stationary Shrödinger equation ĥψ l,q (x) = Eψ l,q (x), after multiplying both sides by e iqx and performing the integral dx, gives the equation for coefficients of decomposition in the Fourier basis 2j + q k