Antonios Pavlidis

Marie-Minerve Louërat

Eric Faehn

Anand Kumar

Haralampos-G Stratigopoulos

Symmetry-based A/M-S BIST (SymBIST): Demonstration on a SAR ADC IP

In this paper, we propose a defect-oriented Built-In Self-Test (BIST) paradigm for analog and mixed-signal (A/M-S) Integrated Circuits (ICs), called symmetry-based BIST (Sym-BIST). SymBIST exploits inherent symmetries into the design to generate invariances that should hold true only in defectfree operation. Violation of any of these invariances points to defect detection. We demonstrate SymBIST on a 65nm 10bit Successive Approximation Register (SAR) Analog-to-Digital Converter (ADC) IP by ST Microelectronics.

I. INTRODUCTION

BIST consists in embedding test instruments into the circuit with the aim to facilitate test. BIST can be functional targeting on-chip measurement of performances or defect-oriented targeting the detection of structural defects. Defect-oriented test has gained importance nowadays as it can help reducing test escapes down to sub-ppm levels. It is also a step towards guaranteeing functional safety if it is capable of detecting latent defects, as well as defects that will be triggered in the context of system operation in the field.

In this paper, we present a defect-oriented BIST strategy for A/M-S ICs, called symmetry-based BIST or SymBIST. SymBIST exploits symmetries into the design so as to build invariant properties whose violation points to defect detection. We demonstrate how SymBIST successfully applies to a 65nm 10-bit SAR ADC IP by ST Microelectronics achieving transparency to the operation, low area overhead, low test time, compatibility with 2-pin digital access mechanisms, and high defect coverage.

Existing works on ADC BIST focus on functional BIST [START_REF] Dufort | On-chip analog signal generation for mixed-signal built-in self-test[END_REF]- [START_REF] Barragan | Practical simulation flow for evaluating analog/mixed-signal test techniques[END_REF]. The main reason for the lack of defect-oriented ADC BIST solutions is the long ADC simulation time, typically in the order of hours, which prohibits a defect simulation campaign. Thanks to the fast test time accomplished by SymBIST, and by using the recent mixed-signal defect simulator Tessent®DefectSim by Mentor®, A Siemens Business [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF], we are able to run very fast in a automated workflow defect simulation for hundreds of defects and compute defect coverage for the entire IP.

The rest of the paper is structured as follows. In Section II, we provide an overview of the SymBIST concept. In Section III, we present the architecture of the SAR ADC IP. In Section IV, we show how SymBIST applies to this IP. In Section V, we specify the defect model and we provide an overview of the defect simulator. In Section VI, we present transient simulations and the defect coverage results. Finally, Section VII concludes the paper. The proposed SymBIST paradigm exploits existing symmetries into an A/M-S design and builds invariant properties that should hold true only in defect-free operation.

The underlying observation is that symmetries are inherent to virtually all A/M-S designs. Inherent symmetries exist thanks to fully-differential (FD) signal processing, complementary signal processing, and replication of identical blocks. Symmetries can also be created artificially with reconfiguration using switches, duplication of blocks, or pseudoduplication of blocks [START_REF] Huertas | Testable switched-capacitor filters[END_REF], [START_REF] Stratigopoulos | Concurrent detection of erroneous responses in linear analog circuits[END_REF]. For example, for node pairs carrying FD or complementary signals we can build an invariance in the form of V 1 + V 2 = α, where V 1 and V 2 are the nodes' voltages and α is a constant, i.e., in the case of FD signals α = 2V cm, where V cm is the common mode voltage [START_REF] Stessman | System-level design for test of fully differential analog circuits[END_REF], [START_REF] Stratigopoulos | An adaptive checker for the fully differential analog code[END_REF]. For identical blocks, duplicated blocks, or pseudoduplicated blocks, we can drive them with the same input and build an invariance in the form of V 1 -V 2 = 0, where V 1 and V 2 are the blocks' outputs.

These invariances can be checked with a window comparator circuit implementing a comparison window [-δ, δ], δ > 0, to account for process, voltage, and temperature variations. If the invariance is violated, i.e. the invariant signal slides outside the window, then this points to defect detection. The parameter δ can be set to k • σ, where σ is the standard deviation of the invariant signal computed by a Monte Carlo analysis and k is set accordingly so as to avoid yield loss.

A high-level abstraction of the proposed SymBIST strategy is illustrated in Fig. 1. The A/M-S IC is divided into purely digital blocks on one side and A/M-S blocks on the other side. We assume that the purely digital blocks are tested with standard digital BIST, i.e. with scan insertion and a combination of stuck-at, bridging, I ddq , and transitional Automatic Test Pattern Generation (ATPG). The A/M-S blocks are divided into three groups. The first two groups include blocks that are FD, perform single-to-FD conversion, provide complementary outputs, appear multiple times, etc. For these blocks invariances exist naturally and the SymBIST strategy applies directly. The third group includes the rest of the blocks. For some of these blocks it may be possible to perform reconfiguration or pseudo-duplication so as to apply SymBIST. We also have the option to perform direct duplication of blocks. For the remaining blocks that are not treated with SymBIST, we need to develop other BIST approaches.

III. CASE STUDY: SAR ADC IP

Our case study is a 65nm 10-bit SAR ADC IP by ST Microelectronics. The top-level architecture is illustrated in Fig. 2. The circuit accepts a FD analog input ∆IN=IN + -IN - and provides a 10-bit digital output D<0:9>. The top-level blocks include: SARCELL: It is the main block of the ADC which implements the SAR algorithm.

SAR Control: It creates 12 pulses P<0:11> used to control the sampling, conversion, and digital output capture phases in the SARCELL.

Bandgap: It creates the required biasing for all ADC blocks. Reference Buffer: It creates the comparison levels VREF<0:32> used during the conversion.

The SARCELL block-level architecture is shown in Fig. 3. It includes the following sub-blocks:

10-bit DAC: The DAC sets the comparison level to which the input is compared at each conversion cycle. It has a resistive plus charge redistribution architecture. As shown in Fig. 4, it is composed of two sub-DACs and a switched-capacitor (SC) array. The two sub-DACs, namely SUBDAC1 and SUBDAC2, have 5-bit digital inputs B<5:9> and B<0:4>, respectively, where B<0:9> is the 10-bit digital output of the ADC during the conversion. The sample-and-hold (S&H) operation that keeps the input signal constant during the conversion is performed within the SC array. SUBDAC1 converts the 5 most significant bits (MSBs) B<5:9>, to comparison levels M+ and M-, while SUBDAC2 converts the 5 least significant bits (LSBs) B<0:4> to comparison levels L+ and L-. The Boolean functions implemented by SUBDAC1 and SUBDAC2 are given by:

M+ = VREF[ 9 i=5 B i • 2 i-5 ] M-= VREF[32 - 9 i=5 B i • 2 i-5 ] L+ = VREF[ 4 i=0 B i • 2 i ] L-= VREF[32 - 4 i=0 B i • 2 i ] (1)
Phase Generator: It controls the timing of the operation by generating the phases for sampling, comparison, conversion, etc.

V cm Generator: It generates the common mode voltage V cm used inside the DAC.

SAR Logic: It controls the conversion process by providing the digital input to the DAC, it stores the result of each comparison, and provides the digital output D<0:9>=B<0:9> once the 10 conversion periods are completed.

Comparator: It compares the two outputs of the DAC and the outcome of the comparison is driven to the SAR Logic block in order to set the corresponding digital bit. It comprises a pre-amplifier, a comparator latch, an RS latch, and an offset compensation circuit for the pre-amplifier.

IV. SymBIST APPLIED TO SAR ADC IP 1) Invariances: Looking into the architecture of the SAR ADC IP, we observe that the two sub-DACs within the DAC are structurally identical, each sub-DAC has complementary outputs, the SC array has symmetrical paths, the pre-amplifier is FD, and the comparator and RS latches have complementary outputs. Based on these observations, we can build the following invariances that hold true for any FD input ∆IN and at every conversion cycle: sgn

M+(i) + M-(i) = VREF[32] L+(i) + L-(i) = VREF[32] (2)
DAC+(i) + DAC-(i) = 2V cm (3)
LIN+(i) + LIN-(i) = 2V cm2 (4)
(Q+(i) -Q-(i)) = sgn (LIN+(i) -LIN-(i)) Q+(i) + Q-(i) = VDD (5)
where V cm2 is the common mode at the outputs of the preamplifier, sgn(•) denotes the sign function, the argument (i) denotes the i-th conversion, and the rest of the signals are annotated in Figs. 3 and4.

2) Test stimulus: The same test stimulus is used for checking all of these invariances. First, the FD input ∆IN stays constant at a DC value which can be set arbitrarily. Second, a 5-bit digital counter is used that generates all possible 2 5 bit combinations at the inputs B<0:4> and B<5:9> of the two sub-DACs. The rationale of this dynamic part of the test stimulus is that all components within the DAC are activated and also the comparator is extensively exercised with various inputs. The components within the bandgap and reference buffer are also activated since during this test all comparison levels VREF[j], j = 0, • • • , 32, are used within the DAC, as shown from Eq. ( 1). The V cm Generator is checked directly with the invariance in Eq. ( 3).

3) Coverage of A/M-S part: With the 6 invariances in Eqs. ( 2)-( 5) and the chosen test stimulus, SymBIST covers all A/M-S blocks of the SAR ADC IP. The purely digital blocks, i.e. SAR control, phase generator, and SAR logic, can be covered with standard digital BIST.

4) SymBIST infrastructure: The SymBIST infrastructure includes the 5-bit digital counter and one window comparator per invariance. Alternatively, we can employ a single comparator and switch it to check invariances sequentially. This choice reduces the area overhead at the expense of test time. A window comparator is connected via switches and buffers so as to be non-intrusive and avoid any performance penalty. SymBIST is embedded without imposing any design modifications in the IP. Moreover, since the test stimulus is digital and the comparator's output is a 1-bit pass or fail decision, SymBIST can be interfaced with a 2-pin digital test access mechanism. Overall, the area overhead of the SymBIST infrastructure is estimated to be less than 5%.

5) Test time: Considering the scenario where the 6 invariances are checked sequentially, the test is completed very fast in 6 * 2 5 • (1/f clk ) = 1.23µsec, where f clk = 156MHz is the clock frequency. This time equals about 16x the time to convert one analog input sample.

V

. DEFECT MODELING AND DEFECT SIMULATOR

We rely on a standard defect model that includes short-and open-circuits across transistor and diode terminals and ±50% variations in passive components, i.e. resistors and capacitors. We use a short defect resistance of 10 ohms. A weak pull-up or pull-down is assigned to each open defect to account for the fact that an ideal open does not exist and, besides, it cannot be handled by a SPICE simulator [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF].

Defects are simulated at transistor-level and in an automated workflow using the Tessent®DefectSim tool by Mentor®, A Siemens Business [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF]. Defects are assigned a relative likelihood of occurrence that is estimated by combining global defect-type likelihoods, i.e. the likelihood of short-circuits is typically higher than the likelihood of open-circuits, and component-specific likelihoods, i.e. the expected component area on the layout, as explained in [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF]. For this reason, we report the L-W defect coverage computed by the tool [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF]. To reduce defect simulation time, we use the stop-on-detection and Likelihood-Weighted Random Sampling (LWRS) options [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF]. When the LWRS option is used, the 95% confidence interval of the L-W defect coverage is also reported [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF].

VI. RESULTS

For our experiment we use a comparison window with δ = 5 • σ, i.e. k = 5, so as to guarantee that yield loss is negligible. Fig. 5 shows the invariance signal in Eq. ( 3) for the defectfree case and for three randomly chosen defect cases within the blocks covered by this invariance, namely the SUBDACs, SC array, and V cm generator. The comparison window is illustrated with the two dashed horizontal lines. The stopon-detection option was disabled for these simulations. The instantaneous glitches are due to the switching operation. A clocked comparator is used to check the invariance with the checks performed when the node voltages are settled, thus no defect detection is flagged when the glitches exceed the range. As it can be seen from Fig. 5, while the defect in V cm generator is detectable during the entire test duration, the defects within the SUBDAC1 and SC array are detectable during specific conversion periods.

Table I shows for the individual blocks of the SAR ADC IP and for its complete A/M-S part the total number of defects, the number of defects simulated, the total defect simulation time, and the L-W defect coverage values achieved using SymBIST including the 95% confidence interval when the LWRS option is used.

As it can be seen, for the entire A/M-S part of the IP, the defect coverage is 86.96% ± 3.67%. As a comparison, for two considerably smaller industrial A/M-S IPs, namely a bandgap and a power-on-reset circuit, the reported defect coverage values are 74% and 51%, respectively [START_REF] Sunter | Using mixed-signal defect simulation to close the loop between design and test[END_REF].

For certain blocks the L-W defect coverage values turn out to be very low. In fact, the absolute defect coverage is much higher, but the undetected defects have a high relative likelihood, thus dominating the L-W defect coverage.

The defect simulation times at first depend on the number of defects simulated, but by using the stop-on-detection option they also depend on the percentage of defects being detected, as well as on the detection time stamps during the test duration. For the entire A/M-S part of the SAR ADC IP, simulating 101 defects took approximately 6060s on a server with 16 cores@1.5 GHz and 128 GB RAM.

Undetected defects should be analysed carefully and it is also interesting to report the percentage of undetected defects that result in at least one specification being violated [START_REF] Gil | Assessing AMS-RF test quality by defect simulation[END_REF]. This is a tedious and time-consuming analysis and is out of the scope of this paper.

VII. CONCLUSIONS

We proposed the defect-oriented SymBIST paradigm for A/M-S circuits, which relies on building signals that by construction are invariant and deviate from the nominal expected value in the presence of defects. SymBIST is demonstrated on a 65nm 10-bit SAR ADC IP by ST Microelectronics achieving a L-W defect coverage of over 85% for the entire A/M-S part. SymBIST presents many advantages, including transparency to the design, very low area overhead, minimum test time, and compatibility with 2-pin digital test access mechanisms.

Fig. 1 :

 1 Fig. 1: High-level abstraction of SymBIST strategy. II. SymBIST PRINCIPLE

Fig. 2 :

 2 Fig. 2: Top-level architecture of SAR ADC.

Fig. 3 :

 3 Fig. 3: SARCELL block-level architecture.

Fig. 4 :

 4 Fig. 4: 10-bit DAC block-level architecture.

Fig. 5 :

 5 Fig. 5: Defect detection by checking invariance in Eq. (3).

TABLE I :

 I L-W defect coverage results with SymBIST.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the financial support in the framework of the ANR EDITSoC project with N o ANR-17-CE24-0014-02. The authors would like also to thank Stephen Sunter from Mentor, A Siemens Business, for the technical advice on the defect simulator.


Aurélien Merkel

Vicent Romero-Garcia

Jean-Philippe Groby

Jensen Li

Johan Christensen

Vicent Romero-García

Unidirectional zero sonic reflection in passive PT -symmetric Willis media

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.

Metamaterials and metasurfaces have emerged as artificial structures that enable one to tailor the wave propagation at will [1]. Whether the field of interest is optics, acoustics, or elasticity, recent efforts have enabled exiting wave phenomena ranging from cloaks of invisibility and unhearability, zero-index behaviour and negative refraction, just to name a few of many counterintuitive effects [2][START_REF] Norris | Proc. R. Soc. A[END_REF][START_REF] Torrent | [END_REF][5][6]. Wave-based diodes and rectifiers enabling one-way sound or light propagation have recently also been put on the map of contemporary metamaterial research [7][8][9][10][11][12].

Waves irradiating onto a passive reciprocal medium usually display symmetric transparency in that the transmittance does not depend on the side at which waves are launched. In contrast, if the elementary building units of the medium involved lack intrinsic inversion symmetry, then the reflectance depends on from which side of the slab waves are irradiated. Commonly, for electromagnetic (EM) waves, such asymmetric response are known to occur in bianisotropic media, which are important in applications comprising unidirectional radiation, single-sided light detection and emission [START_REF] Cheng | Proc. IEEE[END_REF][START_REF] Kong | Proc. IEEE[END_REF][START_REF] Kriegler | [END_REF][16][17].

The intrinsic asymmetry of the involved meta-atoms leads to subwavelength cross-coupling where magnetic states can be excited also by electric fields, and likewise, electric states that can be excited through magnetoelectric coupling. An analogue picture for the case of sound is provided through the cross-coupling between strain and velocity in so-called Willis media [18,19]. Local and non-local coupling in the form of microstructural asymmetry and finite phase changes across inhomogeneities, respectively, define the very nature of Willis coupling in bianisotropic artificial acoustic or mechanical media [20][21][22][23][24][25][26]. To this extent, Willis media have not only emerged to obtain physically meaningful effective parameters, but also for the design of passive acoustic one-way reflectionless systems with the potential to engineer asymmetric steering of sound and to create mechanical structures that suppress echoes from one side, but act reflective from the other.

Asymmetries related to wave propagation have been made also possible in systems disobeying parity P and time-reversal T symmetry, but that are nevertheless symmetric under simultaneous operation. PT symmetry or non-Hermiticity in optics or acoustics, embraces the introduction of loss and gain constituents in new synthetic material such as metamaterials. In optics, gain is an essential ingredient to overcome metamaterial losses but has found increasing importance in one-way waveguides, single-sided scattering, and enhanced sensing at the exceptional point (EP), which marks the transition between an exact and a broken PT -symmetry phase [27]. Specifically, among the phenomena occurring at the EP, unidirectional reflectionless wave propagation, which is also called anisotropic transmission resonance, originates from a non-Hermitian degeneracy of the scattering matrix [28][29][30][31][32]. Interestingly, unidirectional zero reflection has also been found in entirely passive non-Hermitian systems without involving gain, thus posing significant ease for experimental realizations [33][34][35][36][37][38][39][40][41].

In this letter we present a theoretical unification framework that merges the properties of the acoustic constitutive relations containing local and nonlocal Willis coupling with a conventional PT symmetric system representation at the EP in asymmetrically side-loaded waveguides. In details, based on experimental data, we map the one-way reflectionless sound propagation in Willis media onto an ideal, that is, a gain and loss balanced PT symmetric system through a gauge transformation comprising a shift in the P operation associated with an average loss bias. This unconventional strategy to tailor the asymmetric flow of sound further extends the family of non-Hermitian wave physics to synthesized PT matter.

We begin by considering an effective medium containing Willis coupling, hence, the conservation of momentum in an acoustic waveguide is

∇P = - ρ eff ρ 0 c 0 S ∂U ∂t -χ 1 ∂P ∂t , (1)
and the conservation of mass reads

∇ • U = - ρ 0 c 0 S K eff ∂P ∂t -χ 2 • ∂U ∂t , (2)
where P ≡ P /(ρ 0 c 0 ) with the acoustic pressure P , ρ 0 is the mass density and c 0 is the speed of sound in the background medium, U ≡ SU with the acoustic velocity U and S is the cross-section of the cylindrical waveguide, χ 1,2 are the Willis coupling parameters, and ρ eff , K eff represent the effective mass density and bulk modulus, respectively. Because of reciprocity, the relation -χ 1 = χ 2 ≡ χ is always satisfied. The constitutive matrix M eff in one dimension along the z with the time convention e -iωt , has the following form

∂ ∂z Π P U = iω Sρ 0 c 0 /K eff χ -χ ρ eff /(Sρ 0 c 0 ) P U , (3)
with the operator

Π = 0 1 1 0 . ( 4
)
to acquire a form similarity to the constitutive matrix of chiral media for EM waves and Hamiltonian mechanics.

We can find the elements of the constitutive matrix from the scattering coefficients, which can be derived theoretically from the transfer matrix method [42,43] and measured experimenally [44]. The acoustic pressure and velocity located before and after one cell of length L composed of this effective medium can be related via the constitutive matrix M eff as follows

U z=L P z=L = e iωΠM eff L U z=0 P z=0 , (5)
leading to

e iωΠM eff L = 1/(ρ 0 c 0 ) 1/(ρ 0 c 0 ) 1/(SZ w ) -1/(SZ w ) t -r L r R /t r R /t -r L /t 1/t × 1/(ρ 0 c 0 ) 1/(ρ 0 c 0 ) 1/(SZ w ) -1/(SZ w ) -1 , ( 6
)
where t is the complex and reciprocal transmission coefficient, r L,R are the complex reflection coefficients from left L and right R incidence, respectively, and Z w is the acoustic impedance of the waveguide.The effective parameters are then derived from the matrix logarithm of the left-hand side in Eq. (6). In order to realize such effective metamaterial system containing an asymmetric subwavelength cell, we fabricated an acoustic waveguide with two side-loaded resonators tuned at two distinct resonance frequencies as depicted in Fig. 1(a). We take the case of two quarter-wavelength resonators of radius R 1,2 = 5.35 mm having their resonance frequencies at f 1 = 895 Hz and f 2 = 862 Hz, corresponding to the lengths L 1 = 9.141 cm and L 2 = 9.51 cm, respectively and separated by a distance L = 3 cm in the waveguide of radius R w = 1.5 cm. The scattering coefficients are measured in an impedance tube with the four microphone method [44]. Whereas peculiar features do not appear obvious in ρ eff and K eff (in the remainder of the article including figures, we defined ρ eff ≡ ρ eff /(Sρ 0 c 0 ) and 1/K eff ≡ (Sρ 0 c 0 )/K eff ), the eigenvalues ζ 1,2 of the constitutive matrix M eff coalesce at approximately 860 Hz as can be seen in Fig. 1(c) demonstrating the formation of an EP. At this point, we show that the constitutive matrix M eff can be mapped onto a conventional 2×2 PT symmetric system matrix of the form [35,[45][46][47]]

Υ P T = a + i∆γ iκ -iκ a -i∆γ , (7)
where a, ∆γ and κ are real numbers. Similarly, M eff can also be mapped onto an entirely real 2×2 matrix that has the form

Υ R = a + ∆γ κ -κ a -∆γ . ( 8
)
These two PT -symmetric two-level matrices have an EP occurring when ∆γ = ±κ. At the EP of our system, FIG. 2. (color online) Gauge transformation of Eq. ( 9) comprising the elements of the constitutive matrix M eff at the EP. The gauge transformation maps the constitutive matrix onto a PT -symmetric matrix ΥP T + iγ0I, where Re(χ) = 0, marked with vertical gray dotted lines in the gray areas or onto a real matrix ΥR + iγ 0 I, where Im(χ) = 0, marked with vertical gray dashed lines in the blue areas.

Re(K -1 eff ) = Re(ρ eff ), Im(K -1 eff ) = -Im(ρ eff ) and χ is complex as can be seen in Fig. 1(b), hence, at first sight, no form similarity seems to exist between M eff , Υ P T or Υ R . However, the constitutive matrix M eff can be gauged through an unitary transformation [48]

M G = U u M eff U -1 u , (9)
where U u is the propagator matrix

U u = cos(kδ) i sin(kδ) i sin(kδ) cos(kδ) , ( 10
)
in which k is the wavenumber and is real. The matrix U u is unitary, i.e., U -1 u = U † u , ensuring the eigenvalues of any transformed matrix G remain invariant. This unitary gauge transformation maps ideal PT symmetric systems such as Υ P T and Υ R of symmetric P operation, that is, z → -z, which is a mirror operation at the center of the cell onto shifted PT symmetric systems accommodating the generalized P operation z → 2δ -z, which is a mirror operation at the position δ from the center of the cell [49]. We begin by gauging the constitutive matrix M eff through the dimensionless factor kδ in order to obtain fully imaginary Willis coupling, Re(χ) = 0, as highlighted with vertical dotted lines in the gray areas in Fig. 2, which is accompanied by two exact analytical constitutive relationships Re(K -1 eff ) = Re(ρ eff ) and Im(K -1 eff -ρ ef f ) = ±2Im(χ). The resulting gauged constitutive matrix then obtains ±Im(χ) in the off-diagonal terms whereas the real parts of the diagonal terms become Re(K -1 ef f ). Since Υ P T represents an ideal balance of loss and gain, M eff cannot exactly assume the same form. However, in order to accomplish this, we introduce an additional gauge transformation that is biasing the system with an average level of loss γ 0 [33][34][35][36][37][38][39][40][41]. In doing this, the ideal PT -symmetric Hamiltonian is mapped onto an effective but passive Hamiltonian that yields a loss-biased constitutive matrix of the following form

M eff = U u a + i∆γ + iγ 0 iκ -iκ a -i∆γ + iγ 0 U -1 u , = U u (Υ P T + iγ 0 I)]U -1 u , (11)
where I is the identity matrix and the condition for the EP remains as earlier, ∆γ = ±κ. Apart from the above, we can also choose to gauge the system towards fully real Willis coupling such that Im(χ) = 0. In this case, one can see in Fig. 2 marked with vertical gray dashed lines in the blue areas that Im(K -1 eff ) = Im(ρ ef f ) = -iγ 0 and that Re(K -1 eff -ρ eff ) = ±2Re(χ), which correspond to the EP condition for the physical effective system when transformed from Υ R

M eff = U u a + ∆γ + iγ 0 κ -κ a -∆γ + iγ 0 U -1 u , = U u (Υ R + iγ 0 I)]U -1 u . (12)
Thanks to the fact that the signature of PT symmetry breaking continues to exist even when an average loss bias is applied to ideal PT Hamiltonians, the EP embedded in the passive acoustic system will not cease to exist. Hence, our passive system that is described by effective constitutive relations containing the Willis coupling can be mapped onto two different forms of the ideal PT Hamiltonians using an unitary gauge transformation, which conclusively unifies the resulting asymmetry in the acoustic reflections.

In what follows, we show how the EP embedded in the M eff matrix translates into the scattering properties of the sample. The reciprocal scattering matrix that relates the amplitude of the outcoming waves to the amplitude of the incoming waves is here defined as

S = t r R r L t . (13)
Its eigenvalues are ξ 1,2 = t ± (r L r R ) 1/2 and its eigenvectors can be written either as

( √ r R , ± √ r L ) or (± √ r R , √ r L ).
The coalescence of the eigenstates of this scattering matrix is characterized by either r L = 0 or r R = 0 and r L = r R . In Fig. 1(c) we demonstrated that the eigenvalues ζ 1,2 of the constitutive matrix coalesce at 860 Hz and agrees perfectly well with the spectral location of the EP of the associated scattering matrix as it can be seen in Fig. 3(a). The coalescence of the scattering eigenvalues at that point yields reflectionless propagation for left incidence r L → 0 = r R as seen in Fig. 3(b), where experimental data agree very well to theoretical predictions in asymmetrically side-loaded waveguides. We further emphasize that the size of the cell is substantially smaller than the incident wavelength λ, i.e., L < λ/13 and that the unidirectional reflection-less behavior is associated with high levels of absorption within the sample and a consequential low transmittance. For right incidence on the other hand, the wave is mostly reflected with lower level of absorption, denoting a high level of asymmetry in both the absorptions and reflections. Interestingly, the PT conditions for the one-dimensional scattering matrix using the generalized parity operator without considering the loss bias assume now the following form

r L r * R e -4ikδ = r * L r R e 4ikδ = 1 -|t| 2 , (14)
r L t * + r * L te 4ikδ = r R t * + r * R te -4ikδ = 0. ( 15
)
We stress that unidirectional reflectionless sound propagation in principle can take place at other frequencies and at a broad variety of resonator configurations [44].

In summary, we have demonstrated that an EP can appear in the constitutive relations of a passive acoustic Willis media, which can be mapped onto ideal PT Hamiltonians through a gauge transformation and an average loss bias. This EP translates into an unidirectional reflectionless propagation, which is, of primary importance for sound absorption because, if combined with coherent perfect absorption, it results in an unidirectional perfect absorber [50]. Further, our present findings can lead to a deeper insight into the unidirectional invisibility phenomena [28,30], altogether showing how the Willis coupling broadens the possibilities of embracing both worlds of acoustic metamaterials and PT symmetry physics at once to achieve unprecedented control of sound and vibrations.

J.C. acknowledges the support from the European Research Council (ERC) through the Starting Grant No. 714577 PHONOMETA and from the MINECO through a Ramón y Cajal grant (Grant No. RYC-2015-17156). J.-P. Groby and V. Romero García gratefully acknowledge the support from the ANR Project METAUDIBLE No.

FIG. 1 .

 1 FIG. 1. (color online) (a) Picture of the experimental setup with the acoustic waveguide side-loaded by two quarterwavelength resonators of different resonance frequencies. (b) Experimentally derived effective parameters K eff , ρ eff and χ. (c) Theoretical (continuous curves) and experimental (dots) eigenvalues ζ1,2 of the constitutive matrix M eff . The EP where the eigenvalues coalesce is marked with an orange vertical line.

FIG. 3 .

 3 FIG. 3. (color online) (a) Eigenvalues spectrum of the scattering matrix. The EP that originates from the coalescence of the eigenvalues of the scattering matrix, is marked with the crossing of orange lines. (b) Theoretical (continuous curves) and experimental (dots) scattering coefficients of the sample. The EP where rL → 0 is marked with an orange vertical line.

ANR-13-BS09-0003 co-founded by ANR and FRAE and the support from the RFI Le Mans Acoustique (Région Pays de la Loire) PavNat project. This article is based upon work from COST Action DENORMS CA15125, supported by COST (European Cooperation in Science and Technology). JL would like to acknowledge support from Research Grants Council of Hong Kong (Project number 16302218).


Edoardo Sterpetti

Johan Biscaras

Andreas Erb

Abhay Shukla

Crossover to strange metal phase: quantum criticality in one unit cell Bi 2 Sr 2 CaCu 2 O 8+x

published or not. The documents may come

Crossover to strange metal phase: quantum criticality in one unit cell Bi 2 Sr 2 CaCu 2 O 8+x

INTRODUCTION

The determination of the phase diagram of Bi 2 Sr 2 CaCu 2 O 8+x (BSCCO-2212) and in particular the crossover lines between different regions in the vicinity of the superconducting dome is important because of the possible existence of a Quantum Critical Point (QCP) [1][2][3]. This QCP is thought to demarcate a quantum phase transition occurring at 0 K upon variation of a tuning parameter [4][5][6][7][8][9]. The tuning parameter here is the charge carrier density. In BSCCO-2212 the pseudogap to strange metal phase boundary and the strange metal to Fermi liquid crossover, common to the phase diagram of all high temperature superconductors, are not well-established. In the pseudogap phase at low doping, an energy range close to the Fermi level exhibits a paucity of available states and is effectively gapped at low temperatures. The strange metal phase at higher doping shows metallic decrease of resistance with temperature in the normal state but with a linear dependence, even at low temperatures where in conventional metals Fermi-liquid theory predicts a T 2 dependence due to the predominance of electron-electron scattering. The strange metal to Fermi liquid crossover has been rarely studied [10][11][12][13] probably due to the difficulty in fabricating sufficiently overdoped samples and adapted methods. Several methods, spectroscopic or otherwise, have been used to investigate and detect the opening of the pseudogap (occurring at temperature T * ) as a function of doping. Among these we can cite Angle Resolved Photo-Emission Spectroscopy (ARPES) [14], Raman spectroscopy [15], tunneling spectroscopy [16][17][18], c-axis resistivity [19] and Inelastic Magnetic Neutron Scattering (MI-NS) [20].

For a direct identification of the critical regime [8] one needs to measure the crossovers down to very low temperatures approaching 0 K as a function of doping and precisely determine the corresponding sheet carrier density at each doping value. Transport data cannot directly access an eventual QCP masked by the superconducting state unless superconductivity is completely suppressed by a very high magnetic field [21,22]. Indirectly however, the crossover lines found above the superconducting dome can indicate the doping range where criticality could occur. This requires multiple and precise temperature dependent transport measurements around the optimal doping value and as close as possible to the superconducting region as well as an interpretation of the transport data in terms of variations in temperature dependence in different regions of the phase diagram. The position of the QCP is an important point of discussion for theories and experiments seeking to explain this phase diagram and the physics of high temperature superconductors and some even postulate the existence of two such points. Many works [2,8,9] argue that this point is characteristic of the normal state strange metal phase and as such should be found at the extrapolation of the boundaries delimiting this phase in the normal state. Carrier density determination in these compounds is generally based on an empirical and commonly used relation [23] which relates the critical temperature T c to doped holes (p) per Cu atom (holes/Cu). Additionally one can also measure the Hall coefficient R H and determine a Hall carrier density n H = 1/qR H . In strongly correlated electronic systems such as cuprates the Hall estimation does not directly give the value of the carrier density as in conventional metals, but a quantity related to it. It is also temperature dependent confirming that the complex underlying electronic structure cannot be interpreted as a single parabolic band. Nevertheless as we showed earlier [24], 1/qR H at fixed low temperature is a useful relative indicator of carrier density modulation [21,[25][26][27]. In what follows the carrier density is that given by the empirical formula unless specified otherwise.

In a previous work [24] we established a comprehensive phase diagram of two-dimensional (2D) BSCCO-2212 as a function of doping and of disorder by achieving a large variation in doping on each side of the superconducting dome by means of our original space charge doping method [28,29]. In this work we focus on the underlying physics of the intrinsic phase diagram and concentrate on the positioning of the extrapolated quantum critical point near the top of the superconducting dome in the critical region of the phase diagram using a low disorder, 1 unit cell (u.c.) thick sample. It should be noted that other methods have been used to dope surfaces of superconductors by in situ depostion of atoms [30], ionic liquids [31] or ozone annealing [32]. Any electrostatic doping technique is efficient over a thin 2D layer determined by the screening length. In BSCCO-2212, at a thickness of 2.5 u.c. [24], discrepancies can be seen due to decreasing efficiency of the technique over the whole sample thickness, which are avoided in 1 u.c. samples. Transport measurements are affected by disorder which changes both the absolute value of the measured sheet resistance and its dependence on temperature [24]. In general several different samples are used for phase diagram measurements when doping is changed chemically. This involves stoichiometric changes which by definition introduce defects. Space charge doping, as opposed to chemical doping, does not affect the crystal at the atomic level and does not introduce disorder. The defect level in our single sample stays constant over all measurements since no change of stoichiometry is involved. This is an obvious advantage over a changing defect level which can induce changes in scattering rates and mechanisms over and above a change in the physical phase of the sample.

SAMPLE FABRICATION

The superconducting 2D BSCCO-2212 single crystal device is sample E with T c (p opt ) = 80 K. Some measurements from sample D of ref. [24] (maximum critical temperature T c (p opt ) = 81 K) are also included to guarantee consistency. Both devices are 3 nm (1 u.c.) thick and were fabricated on 0.5 mm thick soda-lime glass substrates by the anodic bonding technique [33,34]. This technique, which uses the electrostatic attraction of a surface space charge induced in the glass substrate for the adherence of flakes of layered materials placed on the substrate, allows us to obtain high quality samples with lateral dimensions of about 100 µm. Sample thickness is verified by atomic force microscopy and samples are shaped mechanically with a micromanipulator tip to isolate the area of the uniform required thickness. The large lateral size is essential because we use shadow stencil masks aligned on the sample to thermally evaporate gold contacts in a van der Pauw configuration. We thus avoid chemical degradation and surface contamination from standard electron beam lithography processes, which in our experience degrade few-unit-cell thick BSCCO-2212 crystals, eventually making them insulating [24]. Transport and Hall measurements were performed in a high-vacuum cryostat coupled to an external electromagnet. This set-up allows us to vary the temperature in the 3-400 K range as well as to apply a magnetic field of up to ±2 T. Doping is tuned in the device by means of our original space charge doping method [24,28,29] by monitoring the sheet resistance R S at the doping temperature T d = 350 K while applying a gate voltage V G of the order of 100-280 V between the device and a metallic electrode on the opposite face of the glass substrate. The temperature activates slight ionic mobility of sodium ions inherent to soda-lime glas and the electric field engenders their drift. This creates an ionic space charge at the glass-device interface, consisting of an excess or of a depletion layer of sodium ions. The space charge in turn induces a corresponding mirror charge in the device which can reach very high values, in the region of 10 15 cm -2 because of the intense associated electric field. Lowering the temperature below T d freezes the accumulated space charge at the interface by annulling ionic mobility in the glass substrate, thus fixing the carrier density in the device. The change in carrier density is due to the space charge for the conditions given above and not due to ion intercalation from the glass substrate or eventually oxygen diffusion in the sample [24,28,29]. For a given carrier density the temperature dependence of R S is measured from 330 K down to well below T c . Fig. 1 shows the sheet resistance of sample E as a function of temperature for varying charge carrier density. The curve corresponding to the initial doping level is marked by the star. Negative and positive gate voltages V G were applied to increase and decrease the charge density and therefore the critical temperature T c , allowing us to reversibly [24] sweep between the overdoped regime and the underdoped regime as illustrated by the inset of Fig. 1. For these measurements the doping regime is tuned around the optimally doped region where the T -linear behavior of R S (T ) extends to the lowest temperatures, narrowly separating the pseudogap phase and the Fermi-liquid phase [2,3,35,36].

RESULTS

The strange metal phase is characterized by a linear dependence of R S (T ) [35]. This linear dependence changes in other parts of the phase diagram because of different scattering regimes depending on the relative importance of inelastic and elastic scattering [7]. In BSCCO for example there is a downward deviation in the pseudogap phase from the linear dependance of the strange metal phase. To determine the phase boundaries we proceed as shown in Fig. 2 by detecting the temperature corresponding to deviation from linear dependence of R S (T ) [24]. We use a linear fit for the strange metal region with R S (T ) = R 0 + AT , where R 0 is the residual sheet resistance and A is the slope [2]. We define T * and T m as the temperatures below which the resistance shows a downward (for T * ) and upward (for T m ) deviation from the high temperature T -linear behavior. More precisely, T m in the overdoped region for each doping level was extracted as the temperature above which the R s (T ) measurement could be fitted with a linear component and below which it was fitted with the power law R s (T ) = R 0 + BT m . Though a T 2 behaviour typical of the Fermi liquid is expected in the highly overdoped region, earlier experiments in cuprates have established that in the moderately overdoped phase the upward deviation of R s (T ) below T m also follows such a power law and that the exponent m varies between 1 and 2 depending on the material and the doping level [36][37][38]. In our overdoped devices we observed this m(p) dependence

DISCUSSION

The QCP identified with a quantum phase transition is often situated at the extrapolated junction in the superconducting dome and at zero temperature, of the two phase boundaries delimiting the strange metal region [8,9]. In other words if these phase boundaries can be measured to the lowest possible temperatures close to the superconducting dome, the narrow strange metal region between them will provide an indication of the doping range in which the QCP could exist. By this criterion and following our measurements we identify this critical region to be in the range 0.181-0.195 holes/Cu, inside the dome corresponding to the superconducting region. We can thus conclude that this normal state strange metal region at low temperatures extrapolates to what could be a quantum critical point inside the superconducting dome, in agreement with recent findings [39,40]. We note that this identification has been possible both by the use of high quality low-disorder samples and by accurately measuring the temperature dependence of sheet resistance at a stabilized low rate of 0.3-0.5 K/min. These precautions reduce extrinsic influence which could alter the temperature dependence of R S (T ) and influence the determination of these phase boundaries.

Thus we have been able to distinguish the pseudogap and Fermi liquid related deviations of R S (T ) from linear behavior to below 150 K.

Several spectroscopic measurements [14][15][16][17][18]20] have been carried out on bulk samples with charge carrier density varied by controlling oxygen doping. They have mostly identified the critical region by investigating the underdoped part of the phase diagram. They seem to agree on a linear decrease of T * with doping from the underdoped side to the optimal central part of the dome (p = 0.16 holes/Cu), within the substantial latitude permitted by large error bars. Some of them revealed a subsequent flattening on the slightly overdoped side, followed by an abrupt collapse at larger doping [20] with a QCP estimated at p * = 0.19 holes/Cu. As discussed above, n H = 1/qR H as measured from the Hall effect cannot be directly related to the carrier density in these materials but is a useful relative indicator of carrier density modulation. In contrast with the situation in simple metals R H in this case is temperature dependent. Hall measurements at a few temperatures and doping levels in one of our BSCCO-2212 samples are shown in Fig. 4. The R H vs T curves corresponding to different doping levels in Fig. 4a confirm that the maximum of R H (T) is found at T = 120 K and that the variation of R H with temperature decreases as doping is increased from optimal to the overdoped [26], approaching the Fermi liquid part of the phase diagram. We use R H at 120K to calculate the Hall carrier density n H . Fig. 4b shows the anti-symmetrized transverse resistance R xy as a function of applied magnetic field B at T = 120 K.

The phase diagram of Fig. 3 for samples D and E with carrier density inferred by R H is shown in Fig. 5. The range of critical doping as estimated through 1/qR H is [0.88-1.1]×10 15 cm -2 . Despite the difficulty in relating 1/qR H directly to the doping level, an advantage of this approach is that the error in the determination of 1/qR H is relatively constant (±2-4×10 13 cm -2 ) over the whole doping range, being related only to the determination of the slope of the Hall voltage. In contrast the error in the doping value determined by the empirical formula is maximum in the optimal doping region of the superconducting dome. This is because the dome is flat in this region, so that relatively large variations of the charge carrier concentration produce small variations of critical temperature. This is seen in the size of the horizontal error bars in Fig. 3 and in Fig. 5. Specifically, an error of ± 1 K in the determination of T c results in an error on p = ±0.01 holes/Cu around the optimal doping (0.16 holes/Cu) and of ±0.001 holes/Cu around 0.18 holes/Cu. For higher doping the error becomes negligible. Hence, as the critical region is close to the top of the dome, one could expect a larger error in the determination of the critical doping with the empirical formula. It is then reassuring to see that the crossovers and the critical region as deduced from these two estimations of the doping level are similar.

As mentioned earlier, in BSCCO-2212, ground state electronic structure which is effectively masked in part by the superconducting state, could prove to be complex. In some interpretations the quantum critical point is associated with a Lifshitz quantum phase transition, where as a function of increasing doping the hole-like Fermi surface becomes electron-like at a critical doping of about 0.19 holes/Cu, as observed in ARPES spectroscopy and theoretical calculations [15,41,42].

In another cuprate superconductor, YBa 2 Cu 3 O y , the use of very high magnetic fields to suppress superconductivity has allowed to measure a change in carrier density n H =1/qR H , from n H = 1 + p at high doping to n H = p at low doping [21], as a further confirmation of change of the Fermi surface topology at quantum criticality in cuprate compounds [43][44][45].

Figure 1 :

 1 Figure 1: Space charge doping of sample E. R S (T ) curves of sample E at varying doping levels controlled by the gate voltage V G . The curve corresponding to the initial doping is indicated with a star. The R S curve corresponding to the optimal doping level is green. The device has been doped such that holes are reversibly added and removed. The inset shows the corresponding critical temperatures shown as a function of the doping level as determined using the empirical method. The vertical lines indicate the initial doping value and the extremes attained through space charge doping while the arrows indicate the hole and electron doping ranges achieved.

Figure 2 :

 2 Figure 2: Linear fit of the high temperature part of the sheet resistance R S (T ) for sample E corresponding to the three doping levels indicated at the top left. The green lines are the linear fit. The inset shows the difference between R S (T ) and the linear fit for each curve from which T * and T m can be deduced as the temperatures at which deviation from linearity occurs.

Figure 3 :

 3 Figure 3: Phase diagram of 1 u.c.-thick low disorder BSCCO-2212. Doping dependence of the characteristic temperatures T * (filled symbols),T m (empty symbols) and T c (filled symbols on the superconducting dome) of samples D and E. The dashed vertical line indicates the optimal doping level and the violet band indicates the critical region corresponding to the low temperature strange metal phase.

Figure 4 :

 4 Figure 4: Hall measurements. a, Temperature dependence of the Hall coefficient R H for different doping levels tuned by the external gate voltage V G . b, Transverse anti-symmetrized resistance R xy as a function of magnetic field at T = 120 K. Colors correspond to the same doping level in the two panels.

Figure 5 :

 5 Figure 5: Phase diagram of sample E and D with doping inferred from R H . Characteristic temperatures T * , T m and T c as a function of 1/qR H . The 1/qR H scale is logarithmic to account for an exponential relation between the carrier density determined in this manner and that determined by the empirical formula. The dashed vertical line indicates the optimal doping level, namely the 1/qR H value corresponding to the maximum T c measured. The violet band indicates the critical region corresponding to the low temperature strange metal phase.

CONCLUSIONS

  In this work we identify the critical doping regime in 1 unit cell thick 2D BSCCO-2212 by delimiting the strange metal phase to low temperatures just above the superconducting dome using transport measurements and in-situ electrostatic doping in a high quality low disorder device. Searching for this region in the normal state where the pseudogap, strange metal and Fermi liquid phases appear to meet and by fine-tuning the carrier density in the slightly overdoped regime, where experimental results are scarce, we found it to be centered at about ∼ 0.19 holes/Cu and confined in the range ∼ 0.181 -0.195 holes/Cu. Our transport data allows comparison with spectroscopic measurements and establishes bounds for the range of doping where a quantum critical point could exist though we do not access this region as it is masked by superconductivity. Our results demonstrate the potential of space charge doping to develop fundamental understanding of superconductivity and phase transitions in high-temperature superconductors and open new perspectives in the field of fundamental research.

We thank the Institut des NanoSciences de Paris for access to the electromagnet facility. We acknowledge the Consortium des salles blanches d'Ile de France, M. Rosticher and J. Palomo for access to clean room facilities. This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02, and more specifically within the framework of the Cluster of Excellence MATISSE of Sorbonne Université . * abhay.shukla@upmc.fr


Lamberto Dell'elce
email: lamberto.dell-elce@inria.fr

Daniel J Scheeres

Sensitivity of Optimal Control Problems Arising From Their Hamiltonian Structure

Keywords: Hamiltonian systems, Optimal control, Symplectic structures

First-order necessary conditions for optimality reveal the Hamiltonian nature of optimal control problems. Regardless of the overwhelming awareness of this result, the implications that it entails have not been fully explored. We discuss how the symplectic structure of optimal control constrains the flow of sub-volumes in the phase space. Special emphasis is devoted to dynamics in the neighborhood of optimal trajectories and insight is gained into how errors in the initial states affect terminal conditions. Specifically, we prove that if the optimal trajectory does not satisfy a particular condition, then there exists a set of variations in the initial states yielding a greater error in norm when mapped to the terminal time through the state transition matrix. We relate this result to the sensitivity problem in solving indirect problems for optimal control.

Introduction

The symplectic structure of Hamiltonian systems is responsible for fundamental restrictions that apply to the phase flow [Arnold, 1989, de Gosson, 2011]: the conservation of any full-dimensional volume stated by the Liouville's theorem is arguably the most renown among them, whereas other oftentimes-unintuitive constraints impose limitations on how this volume can be deformed or stretched, e.g., Gromov's non-squeezing theorem, Poincaré-Cartan integral invariants, and the Wirtinger inequality [Gromov, 1985, Maruskin et al., 2009].

Understanding these dynamics provides valuable insight into the propagation of uncertainties in Hamiltonian systems [de Gosson and Luef, 2009], with applications to quantum mechanics and astrodynamics [de Gosson, 2009, Hsiao andScheeres, 2007]. In the context of optimal control, the necessary conditions outlined by the Pontryagin maximum principle reveal the well-known Hamiltonian nature of the problem, where the adjoint variables associated to the enforcement of the equations of movement serve as conjugate momenta [Pontryagin, 1987]. Nonetheless, the implications of this result in the perspective of symplectic structures have not been comprehensively investigated in the literature.

Indirect approaches for optimal control, e.g., shooting method, are of interest because they lend themselves to low-dimensional optimization problems. Nonetheless, their success is often jeopardized by their high sensitivity to variations in the initial states and adjoints, so that direct techniques, e.g., collocation, are used instead. An adequate scaling of the problem can mitigate these effects, but deep understanding of this sensitivity mechanism is missing. The ideal situation for the exploitation of indirect techniques is that any perturbation of the initial conditions is systematically mapped to a 'small' perturbation of the terminal ones.

This work is aimed at assessing how the aforementioned constraints manifest themselves in an optimal control problem and how indirect techniques might benefit from the understanding of this mechanism. Specifically, two questions are tackled: (1) What are the features of symplectic structures when constraints on their terminal states are imposed? (2) What is the sensitivity of optimal trajectories with these constraints with respect to variations in the initial states and adjoints?

After briefly recalling some fundamentals on the Hamiltonial nature of optimal control problems in Section 2, the first question is answered in Section 3, where we focus on Poincaré-Cartan integral invariants. We show that they equal zero when evaluated on the manifold that drives the states to the desired terminal point. The second question is approached in Section 4. Here, we introduce a requirement on the state transition matrix of an optimal trajectory whose fulfillment guarantees that, in the presence of any perturbation of the initial state, the norm of the error on the terminal conditions is smaller than the norm of the perturbation at the initial time. Finally, the outcomes of the paper are illustrated in Section 5 by means of a simple example considering a time optimal planar low-thrust orbital transfer.

2 Hamiltonian nature of optimal control problems Let q ∈ R n and u ∈ U ⊆ R m be a set of dynamical and control variables, respectively. The flow of q is governed by the differential equation q = f (q, u). Consider a maneuver aimed at steering the state vector from some initial conditions, q 0 , to a desired target point, q f , after a prescribed time, t f . The trajectory is required to minimize the integral of the convex cost function g(x, u), so that the optimal control problem is formulated as

q * , u * = arg min x,u t f t 0 g(q, u)dt s.t. : q = f (q, u), u ∈ U ∀ t ∈ [t 0 , t f ], q(t 0 ) = q 0 , q(t f ) = q f . (1)
Denote by H the pre-Hamiltonian, H (q, p, u) = -g(q, u) + p T f (q, u),

where p ∈ R n are the time-varying Lagrange multipliers associated to the constraints q = f (q, u). In the context of optimal control, p(t) is generally referred to as adjoint vector. Without loss of generality, we assume that both q and p are dimensionless in this paper. The Pontryagin maximum principle (PMP) [Pontryagin, 1987] outlines necessary conditions that the solution of Problem (1) is required to satisfy. Specifically, the PMP states that the optimal control, u * , has to maximize the pre-Hamiltonian, i.e., u * (q, p) = arg max u H (q, p, u) ,

and that the trajectory of the state and adjoint variables is governed by the Hamilton's equations

q = ∂ H ∂p , ṗ = - ∂ H ∂q , (2)
where the Hamiltonian, H : R 2n → R, is defined as

H (q, p) = H q, p, u * (q, p) .
Indeed, the solution of Problem (1) has also to satisfy the boundary conditions, q(t 0 ) = q 0 and q(t f ) = q f . Defining the vector of canonical coordinates x = q T , p T T , Eq. ( 2) can be recast in the compact form

ẋ = J ∂ H ∂x , (3)
where

J = 0 n I n -I n 0 n ,
and I n and 0 n denote the n × n identity and null matrices, respectively, The solution of Eq. (3) for given initial conditions, x 0 , is denoted by φ(t, x 0 ). Hence, in the framework of Problem (1), x 0 must be such that

I n 0 n φ(0, x 0 ) = q 0 , I n 0 n φ(t f , x 0 ) = q f .
In this paper we are interested in the flow of the neighborhood of the optimal trajectory. By perturbing the initial conditions with the 2n-dimensional vector δx 0 so that x(0) = x 0 + δx 0 , the resulting trajectory can be expressed as

φ(t, x 0 + δx 0 ) = φ(t, x 0 ) + Φ(t, x 0 )δx 0 + O ||δx 0 || 2 (4)
where

Φ(t, x 0 ) = ∂φ(t, x 0 ) ∂x 0 (5)
denotes the state transition matrix (STM). In the remainder of the paper, the initial condition x 0 in Eq. ( 5) is understood, so that the shorthand notation Φ(t) is used.

The second-order term of Eq. ( 4) is such that

lim ||δx 0 ||→0 O ||δx 0 || 2 ||δx 0 || 2 < ∞,
for a fixed time t. Although this quantity may become large and non-negligible when regarded as a function of time for a fixed δx 0 , it is neglected in this paper. The evolution of the STM is governed by the initial value problem

Φ(t) = J ∂ 2 H ∂x 2 x * (t) Φ(t), Φ(0) = I 2n . (6)
From Eq. ( 6), it is straightforward to verify that d dt Φ T JΦ = 0, so that Eq. ( 6) is an infinitesimal symplectomorfism. For this reason and by noting that

Φ(0) = I 2n×2n ∈ Sp(n) is symplectic, it holds that Φ(t) ∈ Sp(n) ∀ t ∈ [t 0 , t f ],
where Sp(n) = A ∈ R 2n×2n : A T JA = J denotes the group of 2n-dimensional real-valued symplectic matrices.

All results discussed in the remainder of the paper are intimately related to the Hamiltonian nature of optimal control disclosed by Eq. ( 2), which holds true in more general formulations than Problem (1), e.g., non-autonomous system, unspecified maneuvering time, and additional path constraints. However, all developments of this work are built upon Problem (1) to avoid unnecessary notations.

Poincaré-Cartan integral invariants

Integral invariants are quantities that are defined by means of an integral on an arbitrary region of the configuration manifold and that conserve their value for any symplectomorfism of the region itself.

Let Ω(t 0 ) be an arbitrary 2k-dimensional sub-volume in the 2n-dimensional configuration space at time t 0 . The sum of the signed projections of Ω(t 0 ) on all combinations of the symplectic subspaces,

{q i 1 , p i 1 , . . . , q i k , p i k } , 1 ≤ i 1 ≤ . . . ≤ i k ≤ n,
is preserved under any symplectomorfism, which includes the mapping of Ω under the Hamiltonian flow, i.e.,

I k (Ω) = I k (φ(t, Ω)) ∀ t ∈ [t 0 , t f ],
where

I k (Ω) = 1≤i 1 ≤...≤i k ≤n Ω dq i 1 ∧ dp i 1 ∧ . . . ∧ dq i k ∧ dp i k , (7)
Fig. 1 Projection of a surface on the j-th symplectic plane. The dotted region denotes the magnitude of the signed projection. The unsigned projection also includes the dashed region.

and the wedge denotes the exterior product. The quantity I k (Ω) is referred to as Poincaré-Cartan invariant of the volume Ω. An intuitive visualization of Poincaré-Cartan invariants is obtained when a 2-dimensional surface is considered. In this case, the integral becomes

I 1 (Ω) = n i=1 Ω dq i ∧ dp i ,
and it consists of the sum of signed projections of the surface Ω on all symplectic planes (q j , p j ). When subject to the dynamical flow, the surface Ω(t) may fold.

In this case the norm of the unsigned and signed projections can differ from each other, as shown in Figure 1. However, the sum of the signed projections on all symplectic planes (q j , p j ) is preserved. When a 2n-dimensional Ω is considered, the Poincaré-Cartan invariant states that a non-degenerate volume of the phase space is preserved under any symplectomorfism. This result is also referred to as Liouville's theorem, and it can be alternatively deduced by noting that Hamiltonian flows are divergence-less.

Poincaré-Cartan integrals assume an additional compelling feature in the context of optimal feedback control [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF], where for any initial state vector, q 0 , the corresponding adjoint vector is tuned such that q(t f ) = q f . The solution to this problem satisfies the Hamilton-Jacobi-Bellman equation and it defines an n-dimensional manifold on the phase space,

F = x 0 ∈ R 2n : q(t f ) = q f ,
such that the projection of its volume on the space of the dynamical variables q at time t f is equal to zero. As a consequence, all projections of F (and of any arbitrary subset of F) on the symplectic planes are zero at t f , so that the corresponding Poincaré-Cartan invariant is also zero. When mapped back to t 0 , the invariance of the integral imposes a fundamental constraint to the manifold F: any subvolume with non-zero Poincaré-Cartan integral cannot belong to the manifold F. For example, a surface parallel to the j-th symplectic plane cannot reside in F, since all terms of Eq. ( 7) but the j-th one are zero, this giving a non-zero I 1 . The same holds true for any 2k-dimensional volume with its faces parallel to k symplectic planes.

Additional insight is obtained by considering a differential neighborhood of the optimal initial conditions, x * (t 0 ). Denote by Φ qq , Φ qp , Φ pq , and Φ pp the n × n blocks of the STM such that

Φ(t) = Φ qq (t) Φ qp (t) Φ pq (t) Φ pp (t) .
The tangent space of F at point x * (0) is given by

dF = δx 0 : Φ qq (t) Φ qp (t) δx 0 = 0 ,
and it is illustrated in Figure 2. Arbitrary subsets of the n-dimensional space, dF, satisfy the aforementioned requirement of zero Poincaré-Cartan integral. In addition, the first n rows of the STM,

   d T 1 . . . d T n    = Φ qq (t) Φ qp (t) ,
form a complementary basis for dF,

dF ∪ span (d 1 , . . . , d n ) = R 2n .
Any perturbation of the initial conditions with some components in the space spanned by d 1 , . . . , d n yields δq(t f ) = 0.

Symplectic eigenskeleton and contracting directions

Differently of various applications in Hamiltonian mechanics, trajectories of interest in optimal control do not coincide with fixed points of the system, so that the associated STMs are time variant and stability results of linear time invariant systems cannot be exploited to characterize the flow in the neighborhood of the trajectories themselves. For this reason, the developments of this Section are based on the symplectic eigenskeleton introduced in [START_REF] Maruskin | Dynamics of symplectic subvolumes[END_REF] instead.

Denote by Φ(t) the STM associated to the trajectory φ(t, x 0 ). The symplectic eigenskeleton of the transformation φ(t, x 0 ) is defined as the eigenbasis of the STM's Gramian, Φ T (t)Φ(t). The eigenvectors of Φ T (t)Φ(t) occur in pairs, (ξ j , ν j ) such that ν j = Jξ j , j = 1, . . . , n.

In addition, the eigenvalues associated to ξ j and ν j are λ j and λ -1 j , respectively, with 0 ≤ λ j ≤ 1. The vectors {ξ 1 , . . . , ξ n , ν 1 , . . . , ν n } form an orthonormal basis, so that, by defining the matrices Ξ = [ξ 1 , . . . , ξ n ] and N = [ν 1 , . . . , ν n ], it holds

[Ξ N ] Ξ T N T = I 2n . (8)
The symplectic eigenskeleton and the eigenvalues λ j provide insight into how errors in the initial conditions are mapped to final states. For example, consider a small sphere in the phase space encircling x 0 . After applying the Hamiltonian flow, the sphere is squeezed in the directions Φ(t)ξ j and stretched in the directions Φ(t)ν j . The amount of stretch, λ j , is closely related to the finite-time Lyapunov exponents.

The symplectic eigenskeleton, allow to decompose the STM in the form

Φ = R Λ 1/2 0 n 0 n Λ -1/2 Ξ T N T , (9)
where

Λ = diag (λ 1 , . . . λ n ) and R is a symplectic rotation, i.e., R ∈ (Sp(n) ∩ SO(2n)) . (10)
Here, SO(2n) = R ∈ R 2n×2n : R T R = I 2n denotes the special orthogonal group in dimension 2n.

To show that Eq. 10 holds true, we note that the last two factors of Eq. ( 9) are invertible and, by using Eq. ( 8) and the definition of Λ, the matrix R is given by

R = Φ [Ξ N ] Λ -1/2 0 n 0 n Λ 1/2 . ( 11
)
The right-hand term of Eq. ( 11) consists of the product of three symplectic matrices, so that R is also symplectic. In addition, the matrix Φ T Φ can be expressed by means of both Eq. ( 9) and its spectral decomposition, which gives

[Ξ N ] Λ 1/2 0 n 0 n Λ -1/2 R T R Λ 1/2 0 n 0 n Λ -1/2 Ξ T N T = [Ξ N ] Λ 0 n 0 n Λ -1 Ξ T N T .
(12) Pre-and post-multiplying Eq.( 12) by

Λ -1/2 0 n 0 n Λ 1/2 Ξ T N T and [Ξ N ] Λ -1/2 0 n 0 n Λ 1/2 , respectively, yields R T R = I 2n ,
which shows that R is also a 2n-dimensional orthogonal matrix.

As mentioned in the introduction of the paper, indirect techniques would benefit if any perturbation of the initial conditions were systematically mapped to a 'small' perturbation of the terminal states. Because q 0 is given, this requirement can be formulated as

||δq(t f )|| ≤ ||δp(0)|| ∀δp(0) ∈ R n , (13)
The more restrictive notion of omni-directional contraction is now introduced. Specifically, an STM is referred to as an omni-directional contraction if

||δq(t f )|| ≤ ||δx 0 || ∀δx 0 ∈ R 2n ,
where

||δq(t f )|| = ||[I n 0 n ] Φ(t f )δx 0 )|| .
The concept of contraction helps understanding how a favorable formulation of the problem relates to the violating directions defined in Section 3.

Below, we establish a necessary and sufficient condition for an STM to be an omni-directional contraction and, consequently, we provide a sufficient condition for the satisfaction of Eq. ( 13). For this purpose, we recall that any symplectic rotation, R, can be written in the form

R = A B -B A , ( 14
)
where A and B are n × n matrices such that

A T A + B T B = I n , A T B -B T A = 0 n .
Theorem 1 Consider an STM and its decomposition according to Eq. ( 9). Let A and B be the n-dimensional block matrices of the symplectic rotation as defined in Eq. ( 14). Then, the STM is an omni-directional contraction if and only if

C := I n -AΛA T + BΛ -1 B T 0. (15)
In the remainder, C is referred to as contracting matrix.

Proof Equation ( 15) can be proven by expressing the norm of δq(t f ) as

||δq(t f )|| 2 = δq(t f ) T δq(t f ) = δx T 0 Φ T qq (t) Φ T qp (t) [Φ qq (t) Φ qp (t)] δx 0 Hence, the difference ||δx 0 || 2 -||δq(t f )|| 2 is nonnegative if Φ T qq (t) Φ T qp (t) [Φ qq (t) Φ qp (t)] I 2n ⇒ [Φ qq (t) Φ qp (t)] Φ T qq (t) Φ T qp (t) I n (16)
because of the properties of singular values. Substituting Eqs. ( 9) and ( 11) into Eq. ( 16) yields the requirement on the positive definiteness of the contracting matrix defined in Eq. ( 15).

We emphasize that Theorem 1 provides a necessary and sufficient condition, so that if the STM is not an omni-directional contraction, then there exists at least one direction in the phase space, δx 0 , such that its mapping through the STM yields a perturbation δq(t f ) larger in norm than δx 0 , i.e.,

{δx 0 : ||δq(t f )|| > ||δx 0 ||} = ∅.
Hence, if Eq. ( 15) does not hold true, as it is most often the case in real-life applications, the effectiveness of indirect techniques often deteriorates for large t f owing to the generally-exponential increase of the eigenvalues λ -1 j . To gain insight into the concept of omni-directional contraction, consider the case where B = 0 n . In this case, A T A = I n , so that Eq. ( 15) simplifies to I n Λ, which holds true because λ j ≤ 1, j = 1, . . . , n. In this case, all the violating directions d j defined in Section 3 can be written as a linear combination of the contracting directions ξ j . When B = 0 n , it is required that the contribution of the expanding directions ν j is small enough to guarantee the satisfaction of Eq. ( 15).

However, being an omni-directional contraction is neither a property of the system nor of the STM, and the eigenvalues of the contracting matrix are prone to be modified if a canonical change of variables is performed. Block-diagonal transformations are particularly of interest because they do not mix states and adjoints. In addition, it is straightforward to show that these eigenvalues are invariant to symplectic rotations with null off-diagonal blocks, so that the canonical transformation should either include some stretching, i.e., scaling, and/or it does not have to preserve angles. This conjecture is further argued in Section 5.

Example

We consider the extremal flow of a planar time-optimal maneuver in the two-body problem. The Hamiltonian of the problem in equinoctial elements is

H = p 4 µ q 3 1 (1 + q 2 cos q 4 + q 3 sin q 4 ) 2 + F q 3 1 µ 2 i=1 (h i (q) • p) 2
where µ and F denote the gravitational parameter and the maximum thrust-tomass ratio, respectively, and

h 1 (q) =        0 sin q 4 cos q 4 0        and h 2 (q) =               
2 p 1 q 1 1 + q 2 cos q 4 + q 3 sin q 4 cos q 4 + q 2 + cos q 4 1 + q 2 cos q 4 + q 3 sin q 4 sin q 4 + q 3 + sin q 4 1 + q 2 cos q 4 + q 3 sin q 4 0

              
. Table 1 lists the numerical values of the parameters and the initial conditions of the trajectory. Figure 3(a) depicts the eigenvalues of the contracting matrix for the unscaled problem, i.e., the length of the semi-latus rectum is measured in kilometers. The stretching of two violating directions grows up quickly at the very beginning of the maneuver and the condition outlined in Eq. ( 15) is never attained.

Parameter

Gravitational parameter, µ 3.986 • 10 9 km 3 s -2 Thrust-to-mass ratio, F 10 -4 N kg -1 Initial conditions Semi-latus rectum, q 1 10 4 km x-component of the eccentricity vector, q 2 0.75 y-component of the eccentricity vector, q 3 0 True longitude, q 4 0 rad Adjoint vector, p [0, 1, 0, 0] Rescaling the problem consist of applying the canonical transformation

x = s 0 n×n 0 n×n s -1 x,
where s is a diagonal matrix. Figure 3 is nearly an omni-directional contraction over several hours. This suggests that the re-scaled problem exhibit less sensitivity with respect to the non-scaled one.

In addition, the positive definiteness of the contracting matrix is also attained whenever the satellite is close to the apogee of the orbit. Hence, solving a twopoint boundary value problem with a shooting algorithm is supposedly easier if the desired terminal conditions of the maneuver are close to the apogee. This information could be also exploited for choosing when to update the solution during its execution. As an example, Figure 5 depicts the trajectory obtained as a solution of a two-point boundary value problem with terminal conditions equal to q(t f ) = [11 4 km, 0.75, 0, nπ] for an arbitrary odd n (i.e., the terminal state is at apogee). Initial conditions and guess of the initial adjoints are provided in Table 1. The guess of the maneuvering time is set to 50 hours1 . Then, a second set of terminal conditions is chosen on this trajectory such that the final anomaly occurs at perigee (blue squares in Figure 5), and the guess of the maneuvering time is re-scaled accordingly. Hence, the initial adjoint vector solution of the two problems is the same. However, convergence is faster in the first case, as shown in Figure 6. Indeed, a thorough and exhaustive investigation is mandatory to confirm this intuition. We also note that the solver is not able to converge to a feasible solution if states are not scaled.

Conclusion

This paper discussed fundamental constraints on the flow of dynamical systems in optimal control. The symplectic structure of the problem and the imposition of terminal constraints resulted in findings that are not usual in other fields of Hamiltonian mechanics.

Useful resources are readily available by computing the STM of optimal trajectories. First, we identified the local directions orthogonal to the manifold returning the exact terminal conditions. Any perturbation of the initial states with components in this subspace returns a non-zero error in the terminal conditions. Second, after introducing the notion of omni-directional contraction, we derived necessary and sufficient conditions guaranteeing that any perturbation of the initial states is mapped to a 'small' perturbation of the terminal conditions. This condition exposes one mechanism to which an optimal control problem has high sensitivity to errors in the initial adjoint variables.

Fig. 2

 2 Fig. 2 Optimal feedback control surface (shaded region) and violating direction d j .

Fig. 3 Fig. 4

 34 Fig.3Eigenvalues of the contracting matrix as a function of time.

Fig. 5

 5 Fig. 5 Solution of the shooting problem (state trajectory). Blue squares indicate terminal conditions at perigee.

Fig. 6

 6 Fig.6Convergence of the two-point boundary value problems with terminal conditions at apogee and perigee.

Table 1

 1 Parameters and initial conditions for the simulations.

For repeatability purposes, the shooting problem is solved by means of the fsolve function of MATLAB using default settings

Acknowledgements This work was partially supported by the Belgian National Fund for Scientific Research (FNRS) and by the Air Force Office of Scientific Research (AFOSR).

On behalf of all authors, the corresponding author states that there is no conflict of interest.


Paolo Izzo

Aaron J Rosengren

Lamberto Dell'elce

Pini Gurfil

PERIODIC CORRECTIONS IN SECULAR MILANKOVITCH THEORY APPLIED TO PASSIVE DEBRIS REMOVAL

Keywords: Medium-Earth Orbits, Chaos, Resonances, Stability Maps. Post Mission Disposal

Most cartographic stability maps advocated for use in the new passive debris removal ideology based on orbital resonances are obtained through crude averaging methods. This means that from an operational perspective, its not clear where in the osculating space one should actually target to place the satellite on a natural disposal trajectory. It is also not obvious what effects the short-periodic terms may have on these re-entry solutions. We will derive the periodic corrections terms for the dominant perturbations affecting Earth satellites and investigate these considerations.

Introduction

The basic idea in the averaging method is to obtain approximate equations for the system evolution that contain only slowly changing variables by exploiting the presence of a small dimensionless parameter ε that characterizes the size of the perturbation. The tacit assumption is that the perturbing forces are sufficiently weak enough that these approximate secular equations of motion can be used to describe the orbital evolution. The perturbation equations in celestial mechanics, relating the time variation of the orbit parameters to the perturbing accelerations, are nonlinear first-order differential equations of the general form ẋ x x = εg g g(x x x,t),

where g g g(x x x,t) is assumed to be T -periodic in t. Equation 1 is trivially solved when ε = 0, yielding the integrals (Keplerian elements) in the unperturbed problem. The method of averaging consists in replacing Equation 1 by the averaged autonomous system

ẋ x x = ε ḡ g g(x x x), (2)
ḡ g g(x x x) = 1 T T 0 g g g(x x x,t) dt, (3)
in which the average is performed over time, and it is understood that x x x in the integrand is to be regarded as a constant during the averaging process. The basis for this approximation is the 'averaging principle' which states that in the general (non-resonant) case, the short period terms removed by averaging cause only small oscillations, which are superimposed on the long-term drift described by the averaged system.

The removal of time or (what amounts to the same thing) mean anomaly requires computing the quadrature of functions depending implicitly on this variable through the true anomaly. Given a quantity F (α α α, M), defined as a function of the dimensionless time variable M in addition to the other orbital elements represented as α α α, the average is defined by

F (α α α) = 1 2π 2π 0 F (α α α, M) dM, (4)
where the orbital elements α α α are held constant in the integration. Although the average is defined with respect to mean anomaly, it is often more easily calculated by means of the true or eccentric anomaly, using the differential relationships

dM 2π = rdE 2πa = r 2 d f 2πab , (5)
where b = a √ 1e 2 is the semi-minor axis; yielding the equivalent forms for averaging:

F (α α α) = 1 2π 2π 0 F (α α α, M) dM = 1 2πa 2π 0 F (α α α, E)r dE = 1 2πab 2π 0 F (α α α, f )r 2 d f . (6)
Note that r can be expressed in terms of f and E as

r =    a(1-e 2 )
1+e cos f , a(1e cos E).

(7)

The general form of Gauss perturbation equations of the element set H H H (angular momentum vector), e e e(eccentricity vector), and l can be stated as

Ḣ H H = r r r • a a a d , (8a)
ė e e = 1

µ v v v • r r r -H H H • a a a d , (8b)
l = - e µ(1 + √ 1 -e 2 ) H(ê e e • r r r)r r r + (r + p)(ê e e • v v v) θ θ θ - 2 na 2 r r r + ẑ z z • r r r H(H + ẑ z z • H H H) H H H • a a a d + n, (8c)
where r r r and v v v are the position and velocity vectors; θ θ θ = H H H • r r r/H. In this work, using the direct approach of Kozai [START_REF] Kozai | The Motion of a Close garEarth Satellite[END_REF], we will derive the short-period correction terms in the vectorial formulation for the J 2 effect only (contributions due to other perturbing forces will be addressed in future research) and apply these to determine the appropriate initial conditions in the osculating sense for targeting disposal regions depicted in mean element space in stability maps.

Body of paper, give a title

We want to characterize the general form of Equations (8) to include the Earth oblateness disturbing effect. The quadrupole disturbing function arising from an oblate planet can be stated in a general vector expression as follows

R 2 = µJ 2 R 2 2r 3 1 -3(r r r • p p p) 2 , (9)
where J 2 is the dominant oblateness gravity field coefficient, R is the mean equatorial radius of the planet, and p p p is a unit vector aligned with the planet's rotation pole, assumed to be fixed in inertial space. The perturbing acceleration is then given by

a a a 2 = - 3µJ 2 R 2 2r 4 1 -5( p p p • r r r) 2 r r r + 2(r r r • p p p) p p p . ( 10
)
Therefore the perturbation equations in Gaussian form can be stated as

Ḣ H H 2 = r r r • a a a 2 = - 3µJ 2 R 2 r 3 ( p p p • r r r) r r r • p p p, (11)
ė e e 2 = 1

µ v v v • r r r -H H H • a a a 2 = - 3J 2 R 2 2r 4 -1 -5(r r r • p p p) 2 H H H • r r r + 2(r r r • p p p)( v v v • r r r -H H H) • p p p . ( 12
)
The averaged equations of motion have the form

Ḣ H H 2 = 3µJ 2 R 2 2a 3 h 3 ( p p p • ĥ h h) ĥ h h • p p p, (13)
ė e e 2 = - 3nJ 2 R 2 4a 2 h 4 1 -5( p p p • ĥ h h) 2 ĥ h h + 2( p p p • ĥ h h) p p p • e e e. (14)
The differential equations for the periodic term can be written by substracting Eq.( 13) from Eq.( 11)

Ḣ H H sp 2 = Ḣ H H 2 -Ḣ H H 2 = -3µJ 2 R 2 p p p • r r rr r r r 3 + ĥ h h ĥ h h 2a 3 h 3 • p p p, (15)
where sp denotes the short-period perturbation and the orbital elements are taken as constant. Note that the independent variable can be transformed from time to true anomaly by using the relation

dt = r 2 H d f . (16)
The short-period perturbations of the first-order can be obtained as

∆H H H sp = -3µJ 2 R 2 p p p • r r rr r r r 3 + ĥ h h ĥ h h 2a 3 h 3 dt • p p p (17)
where

r r r = cos f ê e e + sin f ê e e ⊥ , (18)
Then the mean value with respect to M can be computed from

∆H H H sp = 1 2π 2π 0 ∆H H H sp dM. ( 19
)
Let us solve the integral in Eq.17

= µ 12H 3 (6 f + 3 sin 2 f + 9e sin f + e sin 3 f ) ê
e eê e e

+ µ 12H 3 (6 f -3 sin 2 f + 3e sin f -e sin 3 f ) ê e e ⊥ ê e e ⊥ (20)
+ µ 12H 3 (-3e cos f -e cos 3 f -3 cos 2 f + 4e + 3) (ê e eê e e ⊥ + ê e e ⊥ ê e e)
and part B B = t 0 ĥ h h ĥ h h 2a 3 h 3 dt = ĥ h h ĥ h h 2a 3 h 3 t = ĥ h h ĥ h h

2na 3 h 3 M (21)
Substituting Eq.( 20) and Eq.( 21) in Eq.( 17) we get

∆H H H sp = -3µJ 2 R 2 p p p • (A + B) • p p p (22)
To obtain the average short period correction, as stated in Eq.( 19), we have to integrate again over the mean anomaly. The core integrals are given by the following equations Eqs.( 23),( 24),(25),(26):

B1 = 1 2π 2π 0 ĥ h h ĥ h h 2na 3 h 3 M = π 2na 3 h 3 ĥ h h ĥ h h (23) A1 = µ 12H 3 1 2π 2π 0 (6 f + 3 sin 2 f + 9e sin f + e sin 3 f ) ê
e eê e edM = 0 (24)

A2 = µ 12H 3 1 2π 2π 0 (6 f -3 sin 2 f + 3e sin f -e sin 3 f ) ê e e ⊥ ê e e ⊥ dM = 0 (25)
Eq.( 24) and Eq.( 25) are equal to zero being the integrals of odd functions over an interval from 0 to 2π

A3 = µ 12H 3 1 2π 2π 0 (-3e cos f -e cos 3 f -3 cos 2 f + 4e + 3) dM = - µ 12H 3 3eX 0,1 0 + 3X 0,2 0 + eX 0,3 0 -4e -3 (26)
Where the Hansen's coefficients of interest are

X 0,m 0 = 1 2π 2π 0 cos(m f )dM = (-e) m (1 + m √ 1 -e 2 ) (1 + √ 1 -e 2 ) m ; m = 1, 2, 3, ... (27)
therefore we finally get the short period correction for the angular momentum vector

∆H H H sp = ∆H H H sp = -3µJ 2 R
The same procedure has to be applied on the eccentricity vector so we substract Eq.( 14) from Eq.( 12) ė e e sp 2 = ė e e 2 -ė e e 2 (30)

The differential equations for the periodic term can be written as

ė e e sp 2 = 3J 2 R 2 4 2 H H H • r r r r 4 - 10 r 4 (r r r • p p p) 2 r r r • H H H - 4 r 4 (r r r • p p p)( v v v • r r r -H H H) • p p p + n a 2 h 4 1 -5( p p p • ĥ h h) 2 ĥ h h • e e e + 2n a 2 h 4 ( p p p • ĥ h h) p p p • e e e (31)
We have to integrate Eq.(31) over time (or true anomaly) to obtain the short-period perturbations of the first-order. To do that, let us rearrange separately some of its terms introducing a new notation that, through the use of dyadic and tryadic expressions, allows to isolate time dependent terms that have to be subsequently integrated.

C1 = r r r r 5 (32) C2 = 1 r 4 (r r r • p p p) 2 r r r • H H H = p • p• rrr r 7 • H H H (33) C3 = 1 r 4 (r r r • p p p)( v v v • r r r -H H H) • p p p = 1 r 4 (r r r • p p p) - µ H 2 (r r re e e + r r rr r r r ) • H H H - µe H (r r r • ê e e ⊥ )U -H H H • p p p = - µ H 2 p • rr r 5 • H H H (e e e • p) - µ H 2 H H H • p• rrr r 6 • p - µe H ( p • rr r 5 • ê e e ⊥ ) p -p • r r 5 ( H H H • p) (34)
In Eq.(34) we considered:

v v v • r r r = r r rv v v -(r r r • v v v)U (35)
where

r r rv v v = - µ H 2 (r r re e e + r r rr r r r ) • H H H (36) r r r • v v v = - µ H 2 r r r • (ẽ e e • H H H) + r r r • ( r r r r • H H H) = µ H 2 [eH(r r r • ê e e ⊥ )] (37)
U is identity dyadic that has the general property: U • p p p = p p p •U = p p p. In Eqs.(33); Eq.( 34) the new notation means that the central part (a dyad or a tryad) is subjected to a scalar or vector product with the surrounding quantities. In this way we can recognize and isolate the core integrals of interest: ∆e e e sp 2 =

I = t 0 r
3J 2 R 2 4 2 H H H • III -10 p • p• I • H H H + 4µ H 2 p • II • H H H (e e e • p) + 4µ H 2 H H H • p• IV • p + 4µe H ( p • II • ê e e ⊥ ) p + 4 ( p • III) ( H H H • p) + M a 2 h 2 1 -5( p p p • ĥ h h) 2 ĥ h h • e e e + 2M a 2 h 2 ( p p p • ĥ h h) p p p • e e e (42)
Then the mean value with respect to the mean anomaly is calculated through the following integrals

I = 1 2π 2π 0 I dM = µ 2 16H 5 -2e 2 X 0,1 0 -e 2 X 0,3 0 - 1 5 e 2 X 0,5 0 -4eX 0,2 0 -eX 0,4 0 -4X 0,1 0 - 4 3 X 0,
+ µ 2 16H 5 -2e 2 X 0,1 0 - 1 3 e 2 X 0,3 0 + 1 5 e 2 X 0,5 0 -4eX 0,2 0 + eX 0,4 0 -12X 0,1 0 - 4 3 X 0,3 0 + 32
III = 1 2π 2π 0 III dM = - µ 2 12H 5 3e 2 X 0,1 0 + e 2 X 0,3 0 + 6eX 0,2 0 + 12X 0,1 0 -4e 2 -6e -12 ê e e ⊥ (44)
= III a ê e e ⊥ (45) The notation we introduced simplify the final expression but, for sake of clarity and to allow an easier implementation, we also report below the expanded algebra written in classical notation summarizing the main steps we did in previous calculations.

IV = 1 2π 2π 0 IV dM = + µ 8H 3 (-eX 0,2 0 - 1 4 eX 0,4 0 -2X 0,1 0 - 2 3 X 0,3 0 + 5 4 e + 8
In order, the first term to consider as it appears in the final expression of Eq.( 47) is:

H H H • III = III a ( H H H • ê e e ⊥ ) (48)
The second term is

1 r 4 (r r r • p p p) 2 r r r • H H H = p • p• rrr r 7 • H H H (49)
after integrating twice, it becomes:

Results

The short period correction curves are obtained correcting only the initial conditions. This corrected set represents the new initial conditions for the singly averaged model. A first-order averaged model, based on the Milankovitch vector formulation of perturbation theory [START_REF] Rosengren | On the Milankovitch Orbital Elements for Perturbed Keplerian Motion[END_REF] which govern the long-term evolution of satellite orbits is summarized in Table 3. where, for completeness, are reported the effects of solar radiation pressure and third-body perturbations, as well as that of the Earth oblateness, studied in this paper. Future work is intended to extend the results here provided to the other perturbations, too.

Table 1. Singly-averaged equations of motion governing solar radiation pressure, planetary oblateness, and third-body gravitational perturbations, where the notation follows from [START_REF] Rosengren | Long-Term Dynamics of High Area-to-Mass Ratio Objects in High-Earth Orbit[END_REF], to which we refer for the omitted details. We test the equations on a reference orbits comparing our result with the following initial orbital elements (semi major axis, inclination, RAAN, argument of perigee and true longitude respectively):

  ê e eê e e + cos f sin f (ê e eê e e ⊥ + ê e e ⊥ ê e e) + sin 2 f ê e e ⊥ ê e e ⊥ (1 + e cos f ) d f

r rr r rr r r r 7 dt = t 0 1 r 4 cos 3 f ê e eê e eê e e + cos 2 f 0 ( 3 f 0 cos 3 8H 3 (2e sin 2 f + 1 4 e sin 4 f + 6 sin f + 2 3 sin 3 f

 02030333 sin f (ê e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + cos f sin 2 f (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê e e) + sin 3 f ê e e ⊥ ê e e ⊥ ê e e ⊥ dt ê e eê e eê e e + cos 2 f sin f (ê e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + cos f sin 2 f (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê e e) + sin 3 f ê e e ⊥ ê e e ⊥ ê e e ⊥ (1 + e cos f ) 5 f + 8e sin 2 f + e sin 4 f + 12 sin f + 4 3 sin 3 f + 12e f ê e eê e eê e e + µ 2 16H 5 -2e 2 cos fe 2 cos 3 f -1 5 e 2 cos 5 f -4e cos 2 fe cos 4 f -4 cos feê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) 5 fe sin 4 f + 4 sin f -4 3 sin 3 f + 4e f (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê 5 f -4e cos 2 f + e cos 4 f -12 cos f -cos f ê e e + sin f ê e e ⊥ )(1 + e cos f ) 2 d f = µ 2 12H 5 (9e 2 sin f + e 2 sin 3 f + 6e sin 2 f + 12 sin f + 12e f )ê e e µ 2 12H 5 (3e 2 cos f + e 2 cos 3 f + 6e cos 2 f + 12 cos f -4e 2 -6e -12)ê e e ⊥ 3 f ê e eê e eê e e + cos 2 f sin f (ê e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + cos f sin 2 f (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê e e) + sin 3 f ê e e ⊥ ê e e ⊥ ê e e ⊥ dt = µ H f ê e eê e eê e e + cos 2 f sin f (ê e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + cos f sin 2 f (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê e e) + sin 3 f ê e e ⊥ ê e e ⊥ ê e e ⊥ (1 + e cos f )d f = µ + 3e f )ê e eê e eê e e e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + e f ) (ê e eê e e ⊥ ê e e ⊥ + ê e e ⊥ ê e eê e e ⊥ + ê e e ⊥ ê e e ⊥ ê e e) );(39);(40);(41) are substituted into Eqs (32);(33);(34) that are subsequently plugged into Eq.(31) to obtain the first-order short-period perturbations.

15 e 2 =

 2 I a (ê e eê e eê e e ⊥ + ê e eê e e ⊥ ê e e + ê e e ⊥ ê e eê e e) + I b (ê e e ⊥ ê e e ⊥ ê e e ⊥ ) = A3(ê e eê e e ⊥ + ê e e ⊥ ê e e)