Sven Dorosz

Thierry Platini

Dragi Karevski
email: karevski@lpm.u-nancy.fr

APS/Fluct-QChains

Keywords:

We study the work fluctuations of two types of finite quantum spin chains under the application of a time-dependent magnetic field in the context of the fluctuation relation and Jarzynski equality. The two types of quantum chains correspond to the integrable Ising quantum chain and the non-integrable XX quantum chain in a longitudinal magnetic field. For several magnetic field protocols, the quantum Crooks and Jarzynski relations are numerically tested and fulfilled. As a more interesting situation, we consider the forcing regime where a periodic magnetic field is applied. In the Ising case we give an exact solution in terms of double confluent Heun functions. We show that the fluctuations of the work performed by the external periodic drift are maximum at a frequency proportional to the amplitude of the field. In the non-integrable case, we show that depending on the field frequency a sharp transition is observed between a Poisson limit work distribution at high frequencies toward a normal work distribution at low frequencies.

I. INTRODUCTION

The study of fluctuations in non-equilibrium small systems has become an active field of research during the last years. The reasons are twofold: on one hand, nanoscaled systems are nowadays quite easily manufactured, opening the door to the emergence of various nano-technologies. On the other hand in the field of non-equilibrium statistical mechanics, where exact results are very few, the discovery of fluctuation symmetries [1,2] expressed by the Gallavotti-Cohen fluctuation theorem [3,4,5,6,7,8] and Jarzynski equality [9] have opened new theoretical perspectives. The fluctuation theorem is a statement on the time-reversal symmetry of the fluctuations of the entropy production along a non-equilibrium path [4]. Whereas this theorem is an asymptotic statement, Crooks derived an interesting identity reading [10,11]

P F (∆S) P B (-∆S) = e ∆S . (1)
∆S is the entropy production for a system driven during a time τ within a forward protocol λ F (t), P F is the distribution of entropy production in the forward process and P B is the distribution associated to the backward process, that is when the system is driven in a timereversed manner. The celebrated Jarzynski equality [9] is easily derived from the Crooks relation. Indeed, utilizing for a thermalised system ∆S = β(W -∆F ), with ∆F = [F (λ F (τ )) -F (λ(0))] and integrating over W we have

e -βW = e -β∆F (2)
which is Jarzynski relation. These relations, initially derived for classical systems, have been extended to the quantal world as well [12,13,14,15,16,17] within various setup and prescriptions for the actual measurement of work performed on a quantum system [18]. As it is clear from the extensivity of the entropy production, the probability to observe a decrease of entropy in a given non-equilibrium trajectory is exponentially small for a macroscopic system. However, for microscopic or mesoscopic systems these untypical trajectories arise with a significant probability and consequently can be observed and measured in actual experimental setup. One may mention in particular experimental tests on the stretching of RNA molecules [19,20], experiments on torsional pendulum [21], on colloidal particles [22,23], on photochromic defect center in diamonds [24]. See Ref. [25] for an excellent review of the experimental related investigations. It is then interesting to calculate or predict the shape of these fluctuations in some explicit models for small system size [26,27]. In this contribution we focus our attention to the work fluctuations of two time-dependent quantum spin chains where the time dependence is due to an externally applied magnetic field. The two models considered are the integrable quantum Ising chain in a time-dependent transverse field and the XX-quantum chain with a longitudinal magnetic field that breaks its integrability. We consider these two cases as archetypical of the integrable and non-integrable situations and we expect that the work distributions will reflect somehow these differences. We check Jarzynski and Crooks relations and compute explicitly the non-equilibrium work distributions associated to several magnetic field protocols. In particular, we study the limiting steady distributions in the important case of periodically driven system.

The paper is organized as follows: in the next section we define the models and protocols for studying work fluctuations. We insist in particular on the definition of work we are using in our study, pointing explicitly as it was done in [18,28] the discrepancies between this definition and the use of a work operator. We present then the numerical results corroborating Jarzynski and Crooks relations before turning to the main part which is the study of work fluctuations within a periodic drift. An exact solution, in terms of double confluent Heun functions is given in the Ising case. The study of the XX-chain is only numerical, but nevertheless evidences are drown out to show clearly the appearance of a sharp transition between an exponential work distribution at high field frequencies toward a gaussian work distribution at lower frequencies. We summarize and discuss our results in the last section.

II. QUANTUM SPIN CHAINS IN A TIME DEPENDENT FIELD

A. Definition of the model and protocol

We study quantum chains with a time-dependent field, in the context of fluctuations of the work performed on them by the time-dependent force, defined through the following Hamiltonian:

H = - 1 2 N -1 i=1 σ x i σ x i+1 - ∆ 2 N -1 i=1 σ z i σ z i+1 - 1 2 h(t) N i=1 σ z i (3)
in particular in the two cases ∆ = 0 and ∆ = 1. The σ's are Pauli's matrices and h(t) is a time-dependend magnetic field applied in the z-direction and leading to the Zeeman term -h(t)M z where M z is the total zcomponent of the magnetization. At ∆ = 0, it corresponds to the Ising quantum chain. In a static transverse field the Ising chain can be mapped after a Jordan-Wigner transformation onto a fermionic problem which can be diagonalised after a suitable canonical (Bogoliubov) transformation. In the thermodynamic limit N → ∞, this model presents a quantum phase transition at h = h c = 1 which is in the 2d classical Ising model universality class. The integrability and non-triviality (in the sense of physical properties) of the model is at the origin of its wide use in many fundamental studies concerning in particular the testing grounds of non-equilibrium quantum statistical mechanics. One may mention for example relaxation properties from an inhomogeneous initial state [29,30,31,32,33,34], entropy production after a local quench [35,36], entropy production in the unique NESS (non-equilibrium steady state) generated from a two temperature initial state [37], or even relaxation properties after a quench through the critical point (basically h is varied from h i > h c to h f < h c ) [38,39].

The case ∆ = 1 corresponds to the isotropic XXquantum chain in a time-dependent magnetic field, which breaks the integrability of the model. With a periodic dependence of the field, this model has to be compared to the quantum kicked Heisenberg chain:

H(t) = - 1 2 L-1 j=0 σ x i σ x i+1 + σ z i σ z i+1 + δ(t)V σ y i σ y i+1 (4)
where the time-dependent perturbation is periodically switched on. This model was studied in the context of the transition from integrability to ergodicity in the thermodynamic limit [40,41,42]. Here instead of the two-body interaction (in the fermionic picture) σ y i σ y i+1 we apply a magnetic field in the z-direction which leads to a manybody fermionic term.

In the following, we consider the distribution of the work performed by the time-dependent field h(t) from the initial time t i to the final time t f . The varying field will lead to transitions between the initial state of the system to a new state at the final time. Here we use as a definition for the work the difference of energies in the final state and initial state,

W if = ∆E if = E f -E i ,
which assumes that we have measurements performed at t i and t f . The work distribution P (W ) is a weighted sum over all initial and final states of the quantum transition probability |(φ f , U φ i )| 2 , where U is the unitary time-evolution operator associated to the time-dependent Hamiltonian H(t), and is given by

P (W ) = i,f δ(∆E if -W )|(φ f , U φ i )| 2 ω φi (5)
where φ i and φ f are respectively initial and final eigenstates of the initial and final Hamiltonian and ω φi = (φ i , ρ 0 φ i ) where ρ 0 ∝ e -βH(0) is the initial (which is supposed to be canonical here) density matrix. It is clear that this definition differs from the introduction of a work operator defined such that its expectation in the state ρ would give the actual performed work:

W = t f ti ds T r{ρ(s) dH ds } (6)
which can be shown to be given by

W = T r{ρ(t i )δH} = T r{ρ(t i )[U + (t f , t i )H(t f )U (t f , t i ) -H(t i )]} (7)
after a short algebra. One is then tempted to define as a work operator δH = H(t f ) -H(t i ), with à the Heisenberg picture of A. However, it is clear that in general the moments (δH) n would differ from the moments W n = dW W n P (W ) with P (W ) defined in (5) and consequently the expectation of the operator e -βδH e -βδH = T r{ρ(t i )e -βδH } (8)

will not necessarily give the Jarzynski result e -β∆F . Nevertheless, one can show that the moments so defined will differ only for n > 2, so that for the average and average square work one can use the work operator δH since then

δH = T r{ρ(t i )δH} = dW W P (W ) (δH) 2 = T r{ρ(t i )(δH) 2 } = dW W 2 P (W ) (9)
with

P (W ) = i,f δ(∆E if -W )|(φ f , U φ i )| 2 ω φi .
The equality of the first moments is trivial. To prove the equality for the second moments, one may notices that for [ρ(t i ), H(t i )] = 0 using the cyclicality property of the trace one has T r{ρ(t i )[ H(t f ), H(t i )]} = 0, which then leads trivially to the equality of the second moments:

dW W 2 P (W ) = i,f (E f -E i ) 2 |(φ f , U φ i )| 2 ω φi = T r{ρ(t i )( H(t f ) -H(t i )) 2 } -T r{ρ(t i )[ H(t f ), H(t i )]} = T r{ρ(t i )( H(t f )-H(t i )) 2 }.
In short, the second moments agree whatever the commutator [ H(t f ), H(t i )] is. At the level of the third moments, n = 3, one may show that a non-vanishing term T r{ρ(t i )[ H(t f ), H(t i )] H(t f )} comes into play leading to a difference between the moments associated to the two preceding definitions of the work.

In order to compute the work distribution one needs to know the unitary time-evolution operator U for the given field protocol, which is a difficult task in general. The numerical study assumes that the true dynamics is well approximated by the step like unitary evolution

U (t, t 0 ) = M-1 n=0 U (t n+1 , t n ) (10)
with

U (t n , t n-1 ) = e iH(tn)∆t ( 11
)
where ∆t is an elementary time increment, meaning that the real function h(t) is approximated by a step function. The true dynamics is obtained by taking the limit ∆t → 0, and M → ∞. In the present study we use typically values of ∆t of order 10 -2 (time is measured in units where = 1 and the spin-spin coupling J x in the xdirection of the spin chains is set to one). The transition probabilities

w i→f = |(φ i , U (t, t 0 ), φ f )| 2
, where φ i (φ f ) are eigenstates of the initial(final) Hamiltonian, are computed numerically using ( 10) and ( 11) on small chains of typically less than 10 spins and with open boundary conditions.

B. Jarzynski and Crooks relations

In order to test our numerical procedure we compare first the free energy differences obtained from the Jarzynski equality (2) in a linear varying field protocol with the equilibrium free energy differences calculated directly from the partition function In figure 1 we have plotted the free energy differences obtained from the Jarzynski equality and from a direct equilibrium calculation as a function of the final field value h at inverse temperature β = 1. The agreement is excellent with a relative deviation which is less than 10 -6 .

F (h(t)) = - 1 β ln T r{e -βH(t) } . ( 12
We have also considered the ratio P F (W )/P B (-W ) = P (W )/P (-W ) for the symmetric protocol h(t) = h o sin(πt/τ ) for t ∈ [0, τ ]. In figure 2 we see explicitly the exponential dependence e βW for h o = 1/2 and 2 for a chain of size N = 8. Again, the agreement of the numerics with the analytical Crooks and Jarzynski relations is fulfilled with a relative deviation which is less than 10 -6 . We have tested also the Crooks relation in the non-symmetrical protocol whose results are presented in figure 3. Again, the exponential dependence e βW is clearly seen on the graph.

C. Stationary distributions in the driven regime

1. Ising quantum chain a. Numerical study We consider now the case where the system is forced with a periodic external field with period τ , maximum value h o and minimum value 0. We are interested in the limiting stationary work distribution (if any) obtained after many periods of the external driving field.

Before going on the driven situation, we concentrate on the fluctuations of the work after one period. In figure 4 we show the distributions obtained on the Ising chain with N = 7 spins for h o = 1/2 and for different periods τ of the forcing field for an initial infinite temperature state β = 0 (equidistribution of the initial eigenstates). In that case, obviously, the average work is vanishing, which is reflected in the symmetry of the work distribution P (W ) = P (-W ). On figure 4 we see clearly that the evolution of the distribution is not monotonic with the forcing period τ . The width of the distribution first increases with τ , reaching a maximum and then decreases with τ . This behavior is clearly seen on the second moment W 2 in figure 5 where we show the variance W 2 β=0 as a function of τ for different maximum field values h o for a N = 7 spin chain. We see that the width of the distribution increases with the amplitude of the perturbation, since then the time-dependent perturbation is more effectively coupled with the spin chain. Moreover, the period τ max associated to the maximum width decreases as the amplitude is increased with, as seen from the numerics, a law τ max ∼ 1/h for large fields enough (h o > 1). Moreover we have seen numerically that τ max is almost independent of the system size. As the field is increased, the maximum value of the variance reaches an asymptotic finite value which is linearly depending on the size N of the chain. One has from the numerics W 2 β=0,max = (N -1)/2. We have also computed the first and second moments of the distribution as a function of τ for different inverse temperatures β. The initial state temperature acts as a scale factor for the moments:

W β = tanh(β/2) δW (13)
where δW , defined by this equation, is the temperature independent part of the average work. As seen on figure 6 one has for the variance

W 2 β -W 2 β ≃ W 2 0 ( 14
)
meaning that the width of the distribution has a very small dependence on the initial state temperature, at least in the considered range of β.

We turn now to the driven situation with many periods. In figure 7 and inverse temperature β = 1. As expected, for large periods the system follows almost adiabatically the variation of the associated equilibrium free energy (the reversible work). Moreover, it is seen as expected from the second principle that the actual average work is always greater than the free energy difference (see the inset).

As the period of the external transverse field is lowered, we see that the deviation from the free energy difference becomes larger and one is no more able to recognize the underlying oscillatory external field. The largest variations and most "chaotic" behavior is obtained at the value τ max previously identified for one period. At higher frequencies the fluctuations of the average work decrease again. b. Exact solution for the driven Ising model With a sinusoidal forcing field h(t) ∝ sin(ωt) it is in fact possible to solve exactly the dynamics of the Ising quantum chain. Indeed, following the lines of Ref. [43], we can put the periodic boundary conditions Hamiltonian (3), after a Jordan Wigner transformation followed by a Bogoliubov one, into a sum of commuting operators:

H = N/2 p=1 Hp ( 15
)
with Hp = (cos φ p -h(t))[a + p a p + a + -p a -p ] ( 16)

-i 1 2 sin φ p [a + p a + -p + a p a -p ] + h(t) (17)
where [ Hp , Hq ] = 0 and where φ p = 2πp/N .

=    h(t) -i sin φ p 0 0 i sin φ p 2 cos φ p -h(t) 0 0 0 0 cos φ p 0 0 0 0 cos φ p    (18)
In the Heisenberg picture, the time evolution matrix U p (t) in the subset p is governed by the equation ( = 1)

i d dt U p (t) = U p (t) Hp (t) (19)
with the boundary condition U p (0) = I. Using ( 18) and (19)

h(t) -i sin φ p i sin φ p 2 cos φ p -h(t) (20)
and U 33 p (t) = U 44 p (t) = e -it cos φp . From the two-coupled first order differential equations, it is easy to obtain the decoupled second order ones. For example, one has for the 11-component

iu ′′ -2 cos φu ′ -(h ′ -i sin 2 φ + ih[2 cos φ -h])u = 0 (21)
with u(0) = 1, u ′ (0) = -ih(0) and similar equations for the other components. Taking for the field the form h(t) = h 1 + h 2 cos(ωt) and putting it in the differential equation ( 21) one gets solutions in terms of double confluent Heun function H D , namely for U 11 p (t) and U 12 p one has

U 11 p (t) = 1 2h 2 (cos φ p -h 1 )e iAt H D (α, β, γ, δ, -iz) -(cos φ p -2h 2 -h 1 )e -iAt H D (-α, β, γ, δ, iz) , ( 22
)
U 12 p (t) = i sin φ p 2h 2 e iAt H D (α, β, -γ, δ, -iz) -e -iAt H D (-α, β, -γ, δ, iz) , (23)
with z = tan(ωt/2),

A t = h2 ω sin(ωt) -t cos φ p , α = γ/2 = 4h 2 /ω, β = 4 ω 2 h 2 2 -|h 2 -h 1 + e iφp | 2 and δ = -4 ω 2 h 2 2 -|h 2 + h 1 + e iφp | 2
, and similar expressions for the other components.

The equilibrium initial state is factorized as ρ(0) ∝ N/2 p=1 e -β Hp(0) and its time-evolution is given by the tensor product

ρ(t) = ρ 1 (t)ρ 2 (t)...ρ N/2 (t) (24)
with

ρ p (t) = U p (t)ρ p (0)U + p (t) . ( 25
)
The average work performed on the system during the time t with the periodic forcing is given by the integral of the magnetization:

W (t) = - t 0 ds h ′ (s)T r {ρ(s)M z } ( 26
)
where

M z = N/2 p=1 M z p , with M z p = a + p a p +a + -p a -p -1
, is the total magnetization operator in the z-direction. Using the factorised form of ρ(t) and the additive structure of M z , it is easy to compute M z (t) as a sum of N/2 independent modes. One obtains

M z (t) = N/2 p=1 tanh βΛp 2 Λ p (h(0) -cos φ p )(2|U 11 p | 2 -1) -2 sin φ p Im U 12 p U 11 p * (27)
with Λ p = sin 2 φ p + (h(0) -cos φ p ) 2 . Together with the solutions of the U p components in terms of Heun functions, this formally solves the problem for the average work. One may notice that the factor tanh( βΛ(0) 2 ) in this expression leads to the previously observed tanh(β/2) in (13) since in that case h(0) = 0 and then Λ p (0) = 1 ∀p and one has for the average magnetization

M z (t) = -tanh β 2 [1 + 2 p cos φ p |U 11 p | 2 + sin φ p Im U 12 p U 11 p * ] ( 28
)
and consequently the average work can be written as in (13): W β = tanh(β/2) δW . The independence of the N/2 modes also explains the linear behavior of the variance observed in the previous section since

W 2 β -W 2 β = N/2 q W 2 q β -W q 2 β ∼ N
, where W q is the work associated to the qth mode. Finally, taking

h 1 = -h 2 = h o /2, such that h(t) = h o sin 2 ( ωt
2 ), we have an exact solution for a situation which closely resembles to the triangular drift treated above numerically. We show the behavior of the average work on figure 8 for different periods τ = 2π/ω. As previously observed, the average work fluctuations are maximum around τ = 3. The inset shows the free energy difference (lower curve) and the average work (top curve) for a large period τ = 50 of the field. In this case the transformation is nearly adiabatic.

XX chain in a longitudinal magnetic field

Finally we present the results obtained in the forced regime for the XX chain. In figure 9 we show the stationary work probability distributions for short and long period τ on a chain of N = 8 spins obtained from an initial infinite temperature state after applying the field over typically several tens of periods (after about 30 periods the work distributions are collapsing nicely on the same curve). We see clearly the appearance of two different regimes depending on the time scale τ . For very short periods τ ≪ 1, the shape of the distribution is very well approximated by an exponential law e -|W |/α /(2α) while for long enough time scale τ its shape is close to a normal law. One may notice on figure 9 the very good collapse of the distributions for different values of the period τ in both short and long periods regimes. We have also seen that this transition survives at finite temperature. The main difference between the infinite and the finite temperature cases being that in the last case the work distribution is no more symmetric but squeezed toward the positive work values as seen in figure 10.

To give evidences of this transition, we have plotted on figure 11 the variance W 2 -W 2 (the average work is vanishing for β = 0) of the distribution P (W ) as a function of the time scale τ . For a system of size N = 8, we observe a sharp increase of the variance from values around 0.25 at time scales smaller than for τ c ≃ 1 to values around 2.5 at larger periods τ where the shape is gaussian-like. Note that the threshold value τ c is depending on the system size, decreasing as the size is increased. On figure 12 we show the field amplitude dependence of the variance W 2 0 of the work distribution for short and long periods of the external drift. One may notice that at long periods, the increase of the variance with the field amplitude is almost linear, while it is more curved at short periods. On figure 13 we see that as the size is increased the variance at large time scales strongly grows with N, while there is almost no dependence on the size at high frequencies. It means that when the period of the external field is to short, the system has no time to follow the perturbation and only transitions between nearby levels are significantly induced. For the density of work w = W/N this leads to a δ(w) distribution in the limit h → 0, since as seen on figure 12, α = W 2 /2 is a vanishing function of the amplitude of the field h, presumably linear. On the contrary, and as the system size is increased, at long periods enough, meaning that the external field drives efficiently the system from its initial level to a new state, transitions between far apart levels are becoming significant. Taking that the variance seems to be const.hN 2 from figures 12 and 13, the distribution π(w) in that case behaves as e -const.w 2 /h . Finally, we have analyzed the temperature dependence of the work distribution P (W ) in the large τ regime. As seen on figure 14, the average work after a first increase with β seems to finally saturate at a value close to 0.44. We have checked that this shape is compatible with a tanh(µβ) behavior. Contrary to the almost temperature independent variance of the Ising case distributions, we see here that the variance of the distribution is decreasing(increasing) with β(temperature), which is what one would normally expect.

III. SUMMARY AND DISCUSSION

In this study we have presented the numerical and analytical results we have obtained for the work distribution of small quantum systems driven by an ex- ternal field. We have considered two different models, namely the integrable Ising quantum chain in a timedependent transverse field and the XX quantum chain in a time-dependent longitudinal magnetic field. In this last model, the presence of the longitudinal field breaks the free fermionic structure of the chain, while the transverse field in the Ising case preserves the free particles structure which obviously leads to integrability. In a first stage, we have checked the validity of both quantum Jarzynski and quantum Crooks relations that were largely discussed in section 2. After this initial check, validating our numerical approach, we focused our attention to the periodic driven situation. In the Ising chain, exploiting its free fermionic structure, we solved exactly the unitary dynamics of the system in terms of somehow complicated double confluent Heun functions. The main features observed are, as expected, that for a slow enough process, that is large period of the oscillations of the field compared to the coupling constant, we recover the adiabatic situation where the work slightly differs, from above, from the equilibrium free energy difference. This is a statement of the second principle. As the process is fastened, the average work starts to deviate significantly from the free energy difference with fluctuations in time that are growing as the frequency of the field is increased. Nevertheless, at very high frequencies, the work fluctuations decrease again to zero. In this case, the variation of the field is so fast that the system is not able to follow it anymore. The maximum amplitude of the work fluctuations is obtained at an intermediate period which is of order of the inverse field amplitude, at least for large field enough. This threshold is understood as a dynamical resonance. In the Ising chain, one is not able to observe a steady work distribution. The behavior of the XX chain in a longitudinal periodic field is quite different. We observe numerically, that the longtime work distribution reaches a steady shape which is strongly dependent on the frequency of the field. Indeed, at high frequencies the distribution is well fitted by an exponential curve, e -|W |/α /(2α), where α ∼ h is almost size-independent. Lowering the frequency, there is a sharp transition toward a Gaussian like behavior for the work distribution with a variance which seems from the numerics to be proportional to hN 2 , . In order to understand the transition between these two limit distribution, one has to realize that the energy band of the chain has a width δ E of order N (since the typical coupling J is set to one). The periodic perturbation introduces a typical energy scale ω ≡ 2π/τ . For ω > δ E , there is no resonant coupling of the system with the periodic forcing. But as soon as ω < δ E resonant coupling leads to the appearance of resonant peaks, that is sharp increases of the transition probability in the regions W = ±ω and integer multiples of ω. The typical width of those peaks is of the order of the amplitude h of the perturbation. Consequently, the deviation to the initial exponential distribution starts at periods τ > O(2π/N ). As the period is increased further, the first peaks moves toward the center of the distribution and new resonant peaks enter into play from the boundaries. Finally, the superposition of these resonant peaks for τ > O(1/h) leads to a new limit distribution of the work. So the transition region between the two limit distributions is O(2π/N ) < τ < O(1/h). Finally, one may remark that the change of shape from exponentiallike to Gaussian-like could be linked to the integrability of the model. Indeed, at very short periods, the system is not able to follow the external perturbation and its behavior is governed by the initial integrable Hamiltonian (XX-chain without field) while at larger periods the system feels effectively its non-integrability, possibly leading to a Gaussian distribution of level spacing (which is our work). To confirm eventually this scenario one needs to

  FIG.1: Free energy difference obtained from the Jarzynski equality compared to the exact value for the Ising chain with N = 7 as a function of the terminal field value at β = 1 and h(t) = t.

FIG. 2 :FIG. 3 :

 23 FIG. 2: Crooks fluctuation ratio in the symmetrical case for the XX chain at ∆ = 1 with N = 8 at β = 1 and τ = 1.

FIG. 4 :

 4 FIG. 4: Work probability distribution for the Ising model at with N = 7 spins at β = 0 for ho = 1/2.

5 FIG. 5 :

 55 FIG. 5: W 2 for the Ising model with N = 7 spins at β = 0 as a function of the period τ .

FIG. 6 :

 6 FIG. 6: Variance of the work distribution for the Ising model with N = 7 spins at ho = 1 as a function of the period τ for β = 0, β = 1/2 and β = 1.

FIG. 7 :

 7 FIG.7: Time evolution of the average work for the Ising chain with a triangular periodic transverse field at β = 1 and maximal amplitude ho = 1/2 with periods τ = 1 (circles), τ = 3 (squares) and τ = 5 (triangles). The inset shows the free energy difference (lower curve) and the average work (top curve) for a large period τ = 50 of the field. In this case the transformation is nearly adiabatic.

FIG. 8 :

 8 FIG.8: Time evolution of the average work for the Ising chain with a sine square periodic transverse field at β = 1 and maximal amplitude ho = 1/2 with periods τ = 1 (circles), τ = 3 (squares) and τ = 4 (triangles). The inset shows the free energy difference (lower curve) and the average work (top curve) for a large period τ = 50 of the field. In this case the transformation is nearly adiabatic.

FIG. 9 :

 9 FIG. 9: Work probability distribution for the XX quantum chain with N = 8 spins at β = 0 for ho = 1/2.

FIG. 10 :

 10 FIG. 10: Work probability distribution for the XX quantum chain with N = 8 spins at β = 1/2 for ho = 1/2.

FIG. 11 :FIG. 12 :

 1112 FIG. 11: Variance of the work distribution for the XX quantum chain with N = 8 spins at β = 0 and β = 1 for ho = 1/2 as a function of the field period τ .

2 >FIG. 13 :

 213 FIG.13: Second moment of the work distribution for the XX quantum chain with ho = 1/2, at β = 0 as a function of the system size N for τ = 0.1 (squares) and τ = 10 (circles).

FIG. 14 :

 14 FIG.14: First (circles), second (squares) and variance (triangles) of the work distribution for the XX quantum chain with N = 8 spins at τ = 10 at ho = 1/2 as a function of the inverse temperature β.

push further this investigation.


Socio-economic evaluation of UK agri-environmental policy: lessons trom an ESA case study

Sarah Skerratt

To cite this version:

  Sarah Skerratt. Socio-economic evaluation of UK agri-environmental policy: lessons trom an ESA case study. Gestion des exploitations et des ressources rurales : Entreprendre, négocier, évaluer, 31, INRA, 437 p., 1998, Etudes et Recherches sur les Systèmes Agraires et le Développement, 2-7380-0843-7. hal-02835348


Thomas Mietton

Luca Rizzi

BRANCHING GEODESICS IN SUB-RIEMANNIAN GEOMETRY

Keywords:

come

Introduction

A metric space is branching if there exist minimal geodesics starting from the same point which follow the same path for some initial time interval, and then become disjoint. Common examples are found in Finsler geometry (e.g. R 2 with the sup norm), or on graphs. On the other hand, it is well-known that Riemannian manifolds and Alexandrov spaces with curvature bounded from below are non-branching.

We are interested here in sub-Riemannian spaces, a large class of metric structures generalizing Riemannian geometry where a metric is defined only on a subset of tangent directions (cf. Section 2 for precise definitions). Several questions concerning geodesics, which are trivial in Riemannian geometry, become hard open problems in the sub-Riemannian setting. For example it has been only recently proven in [START_REF] Hakavuori | Non-minimality of corners in subriemannian geometry[END_REF] that sub-Riemannian geodesics cannot have corners, but it is not yet known whether geodesics are C 1 , see for example [START_REF]Singulières minimisantes en géométrie sous-Riemannienne[END_REF].

To provide further motivation for our contribution, let us mention that there is an on-going effort in trying to define a suitable concept of lower curvature bound for sub-Riemannian spaces, in particular in relation with the synthetic approach à la Lott-Sturm-Villani (cf. for example [BR19, BKS18, BKS19, BG17, Mil19] and [Vil19, p. 1127-43]). Since the existence of branching geodesics causes difficulties in the study of optimal transport and of spaces satisfying synthetic curvature bounds, it is important for further progress in the theory to understand whether sub-Riemannian structures can exhibit such a phenomenon.

In this paper we show that sub-Riemannian (normal) geodesics can branch, adding this phenomenon to the list of remarkable features of sub-Riemannian geometry. Even though normal geodesics are obtained through the action of a Hamiltonian flow, and are therefore smooth, they are not uniquely characterized by their jet at some point. To our best knowledge, it is the first time that this fact is observed.

1.1. Branching and magnetic fields. We describe succinctly the branching phenomenon through an example, of the same basic type as Montgomery's original construction of abnormal minimizers. We will exploit the connection between the equations of motion of a particle in a magnetic field and sub-Riemannian geodesics, pointed out in [START_REF] Montgomery | Isoholonomic problems and some applications[END_REF][START_REF] Montgomery | Abnormal minimizers[END_REF].

Take the distribution on R 3 defined by the kernel of a one-form ω = dz-A(x, y)dy, and consider the sub-Riemannian metric given by the restriction of dx 2 + dy 2 . Let B = ∂ x A be the magnetic field associated with ω, that is ω ∧ dω = -B(x, y)dx ∧ dy ∧ dz. Notice that, without changing B, we can alter A in such a way that A(0, y) = 0 (gauge freedom), and thus the straight line γ 0 (t) = (0, t, 0) for t ∈ R is horizontal.

It is well-known that abnormal paths are precisely the horizontal curves contained in the zero-locus of B, while normal geodesics are those whose projection (x(t), y(t)) on the xy-plane satisfies (1)

κ(t) = λB(x(t), y(t)),
for some λ ∈ R, and where κ(t) is the curvature of (x(t), y(t)), which we assume to be parametrized with constant speed. In particular, the ODE corresponding to (1) describes the motion of a particle with charge λ under the action of the magnetic field B(x, y) normal to the plane. Choose a smooth potential A(x, y) such that B(x, y) = x for y < 0 and B(0, y) > 0 when y > 0. The zero-locus of B coincides in this case with x = 0 when y < 0. In particular, γ 0 (t) = (0, t, 0), for t < 0, is an abnormal geodesic. Since B(0, y) = 0 for y < 0, the curve γ 0 (t) for t < 0 satisfies also (1) for any λ ∈ R, so it is also normal. We can now extend such a curve to a normal geodesic γ λ (t) for t ∈ R, by solving (1) for different values of λ ∈ R. Of course, γ 0 corresponds to the straight line but, from the fact that B(0, y) > 0 for y > 0, the curve γ λ must have nonvanishing curvature for small non-zero λ, hence a branching phenomenon occurs at the origin. Moreover, from what we will show in Proposition 11, the projection on the xy-plane of the trajectories γ λ contains an open neighborhood of the positive yaxis and those trajectories are all distinct for λ sufficiently small. From the physical viewpoint, this phenomenon corresponds to particles having different charges, which "spray out" following different trajectories under the influence of the magnetic field.

In Section 4, we show an explicit construction of such A(x, y), obtained by gluing a flat Martinet structure with the standard Heisenberg one.

Strictly abnormal branching.

In this note we only consider the branching of normal geodesics, and we do not cover the possible branching of strictly abnormal ones (cf. Section 2). It is easy to produce sub-Riemannian structures with branching abnormal paths. For example, consider a degenerate Martinet-type structure in a three-dimensional space, whose Martinet surface itself branches. Such a structure cannot verify the usual non-degeneracy condition, cf. [Mon02, Sec. 3.2]. The Liu-Sussmann local minimality result for abnormal paths does not apply [START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF], and we are not able to prove that these paths are geodesic (i.e. length-minimizing curves).

Open problem. Find an example of branching strictly abnormal geodesics.

Structure of the paper. To make the paper self-contained, in Section 2 we recall some basic facts in sub-Riemannian geometry, following [START_REF] Agrachev | A comprehensive introduction to sub-Riemannian geometry[END_REF]. Sub-Riemannian branching is then discussed in Section 3, and put in relation with the well-known normal/abnormal duality of sub-Riemannian geodesics. The explicit construction of an example, as described in the abstract, is done in Section 4. We conclude by building the most general case of multiply-branching normal geodesics, in Section 5.

Sub-Riemannian geometry

A sub-Riemannian structure on a smooth n-dimensional manifold M , where n ≥ 2, is defined by a set of m global smooth vector fields X 1 , . . . , X m , called a generating frame. The distribution is the possibly rank-varying family of subspaces of the tangent spaces spanned by the vector fields at each point

D x = span{X 1 (x), . . . , X m (x)} ⊆ T x M, ∀ x ∈ M.
The generating frame induces an inner product g x on D x such that

g x (v, v) := min m i=1 u 2 i | m i=1 u i X i (x) = v , ∀v ∈ D x .
We assume that the structure is bracket-generating, i.e., the tangent space T x M is spanned by the vector fields X 1 , . . . , X m and their iterated Lie brackets at x. A horizontal curve γ : [0, 1] → M is an absolutely continuous path such that there exists a

control u ∈ L 2 ([0, 1], R m ) satisfying (2) γ(t) = m i=1 u i (t)X i (γ(t)), a.e. t ∈ [0, 1].
This implies that γ(t) ∈ D γ(t) for almost every t. Notice that the control u in (2) is not unique, but one can always find a unique minimal control, i.e. the one such that g x ( γ(t), γ(t)) = |u(t)| 2 for a.e. t ∈ [0, 1]. We define the length of γ as

(γ) = 1 0 g( γ(t), γ(t))dt.
The sub-Riemannian (or Carnot-Carathéodory) distance is defined by:

d SR (x, y) = inf{ (γ) | γ(0) = x, γ(1) = y, γ horizontal}.
By the Chow-Rashevskii theorem, under the bracket-generating assumption, between any two points x, y ∈ M , there exists a horizontal path, and therefore the sub-Riemannian distance is well-defined. Furthermore, one can prove that it is continuous and its induced topology is the same as the manifold topology. We remark that this definition of sub-Riemannian metric, based on the concept of global generating frame, includes the classical constant-rank case, see [ABB19, Section 3.1.4].

In place of the length, it is convenient to consider the energy

J(γ) = 1 2 1 0 g( γ(t), γ(t))dt.
Horizontal trajectories minimizing the energy with fixed endpoints are exactly paths that minimize the length, parametrized with constant speed. We call a minimizing geodesic between two points x and y in M a horizontal path γ : [0, 1] → M , with γ(0) = x and γ(1) = y, that minimizes the energy among all horizontal paths sharing the same extremities. The term geodesic, instead, denotes the more general class of horizontal paths that are minimizing geodesics locally around each of its points.

2.1. Characterization of geodesics. Let x ∈ M . By Cauchy-Lipschitz, there exists an open set U ⊆ L 2 ([0, 1], R m ) such that, for all u ∈ U, the Cauchy problem:

γu (t) = m i=1 u i (t)X i (γ u (t)), γ u (0) = x admits a solution γ u defined on [0, 1]. The end-point map E x : U → M is then E x (u) = γ u (1)
. The end-point map is weakly continuous and differentiable. The problem of finding a minimizing geodesic between two points x and y is then the problem of minimizing the functional J (seen as a smooth functional defined on U) under the constraint E x (u) = y. By the Lagrange's multipliers rule, if γ is a minimizing geodesic between x and y, and u is its minimal control, then there exists

λ 1 ∈ T * y M and ν ∈ {0, 1}, with (λ 1 , ν) = 0, such that (3) λ 1 • D u E x = ν u, • ,
where • denotes the composition of linear maps and D the (Fréchet) differential. Any path whose minimal control verifies (3) with ν = 0 is called abnormal, or singular.

A path verifying (3) for its minimal control with ν = 1 is called normal. Notice that the case ν = 0 means that u is a critical point of the end-point map.

The covector λ 1 can be interpolated for times t ∈ [0, 1] yielding a lift of the curve γ in the cotangent bundle. In particular in the normal case, the lift λ : [0, 1] → T * M solves a Hamiltonian differential equation. To state this fact more precisely, let us first define some objects: if X is a vector field on M , we can associate to it a function h X : T * M → R defined by h X (λ) = λ, X . In turn, if h is a function on the cotangent bundle, its associated vector field h is defined by σ(•, h) = dh, where σ is the canonical symplectic form on the cotangent bundle, which can be expressed in coordinates as σ = n i=1 dp i ∧ dq i . Finally, for a sub-Riemannian structure in M given by a generating family as above, we define the Hamiltonian H : T * M → R by

H(λ) = 1 2 m i=1 h X i (λ) 2 = 1 2 m i=1 λ, X i 2 .
The following result is an immediate consequence of the characterization of energy minimizers by the Lagrange multipliers rule, or can also be seen as a version of the Pontryagin maximum principle in this setting, cf. [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF].

Theorem 1. Let γ be a horizontal path minimizing the energy between x et y, and let u be its minimal control. Then there exists a Lipschitz path λ :

[0, 1] → T * M lifting γ -that is for all t ∈ [0, 1], λ(t) ∈ T * γ(t) M -such that λ(t) = m i=1 u i (t) h X i (λ(t)) a.e. t ∈ [0, 1].
Moreover, one of those two conditions is verified:

(N) for all t ∈ [0, 1] and all i = 1, . . . , m, it holds u i (t) = h X i (λ(t)), that is λ is solution of the differential equation λ(t) = H(λ(t)); (A) for all t ∈ [0, 1] and all i = 1, . . . , m, it holds h X i (λ(t)) = 0 and λ(t) = 0.
The conditions (N) and (A) correspond to the normal and abnormal cases of Lagrange multipliers rule, with λ(1) corresponding to the multiplier λ 1 in (3).

Remark 2. From (3) the set of normal Lagrange multipliers of a path is an affine space over the linear space generated by its abnormal ones. The same property holds for the corresponding lifts.

The Hamiltonian characterization in the normal case allows us to define an exponential map exp x : T *

x M → M , where exp x (λ) = π • e H (λ), where t → e t H denotes the one-parameter group of diffeomorphisms on the cotangent bundle given by the Hamiltonian flow. In other words exp x (λ) is the extremity at time 1 of the normal geodesic whose lift verifies λ(0) = λ, that is parametrized by constant speed equal to 2H(λ). Normal paths are locally length minimizing, and hence are geodesics. We assume that (M, d SR ) is complete, so that H is a complete vector field.

Note that if a path is normal (resp. abnormal), any smaller segment is also normal (resp. abnormal) as the restriction of the lift verifies the same conditions.

The lift of a minimizing path given by Theorem 1 is not necessarily unique and therefore the same horizontal path can be normal and abnormal at the same time. We will call a path strictly normal if it is normal and it does not admit an abnormal lift, and strictly abnormal if it is abnormal and it does not admit a normal lift.

As a final remark, if a path is strictly normal, then its normal lift is unique. Indeed, if λ and µ were two distinct normal lifts, then λ(1) -µ(1) would be an abnormal multiplier for this path (cf. Remark 2).

Branching geodesics

A natural question is whether strictly normal paths can contain non-trivial abnormal subsegments. We will first show that this behaviour is linked to the occurrence of branching normal geodesics, and moreover that such an occurrence is actually equivalent to a jump in the rank of the differential of the end-point map. Then, in the next section, we will show a simple and natural example of this phenomenon.

First let us define precisely what we will call branching here.

Definition 3. A normal geodesic γ is branching at time t ∈ (0, 1) if there exists a normal geodesic γ such that γ| [0,t] = γ | [0,t] and γ| [0,t+ε] = γ | [0,t+ε] for all ε > 0.

Let γ be a normal geodesic. For t ∈ [0, 1] we define the set

Π t = {λ ∈ T * γ(0) M | γ(s) = π • e s H (λ) ∀s ≤ t}.
This set is a non-empty affine space corresponding to the initial normal covectors of the path γ| [0,t] given by Theorem 1. The set Π t is an isomorphic image of the set of normal Lagrange multipliers of γ| [0,t] and, from Remark 2, its dimension is the corank of the path γ| [0,t] , defined as the corank of the application D ut E x , where u t is the minimal control of γ| [0,t] . The function t → Π t for t ∈ [0, 1] is nonincreasing for the inclusion order and it is thus piecewise constant with some possible jumps where its dimension decreases. Definition 4. The corank function of γ is the function that associates to a time t ∈ [0, 1] the corank of γ| [0,t] , that is the function t → dim Π t . We say that γ is rank-jumping (or corank-jumping) at time t if there is a discontinuity in the corank function for this time.

The corank function is nonincreasing and piecewise constant and moreover, from the lower semicontinuity of the rank and the C 1 regularity of the end-point map, it is left-continuous. We say that γ is rank-jumping (or corank-jumping) at time t if there is a discontinuity in the corank function for this time.

We can now state our theorem linking the phenomenon of branching normal geodesics with rank jumps, which at this point is very elementary to prove.

Theorem 5. A normal geodesic γ branches at time t ∈ (0, 1) if and only if it is rank-jumping at time t.

Proof. Assume γ branches at time t. Then the branching geodesic γ has an initial covector λ such that λ ∈ Π t but λ / ∈ Π t+ε for all ε > 0, which means the rank jumps at t. Conversely, if the rank jumps at time t, there is a covector λ in Π t that is not contained in Π s for s > t, and the path γ defined by γ (t) = π • e t H (λ ) branches with γ at time t.

An immediate consequence is that a normal path can only branch a finite amount of times (up to the maximal corank of a path, which is the corank of the distribution), and furthermore γ| [0,t] must be abnormal.

If γ is strictly normal, its corank is 0 and we have the following corollary, corresponding to the situation encountered in the example from next section. Corollary 6. A strictly normal geodesic γ is branching for some time t ∈ (0, 1) if and only if it contains a non-trivial abnormal subsegment that starts at time 0. In particular if t is the last branching time, γ| [0,t] is a maximal abnormal subsegment.

In this situation, if γ branches at time t ∈ (0, 1), then it branches in a whole family of distinct normal paths, parametrized by the abnormal Lagrange multipliers of the abnormal subsegment. To be more precise, let A ⊂ T * γ(t) M be the set of abnormal Lagrange multipliers associated with the maximal abnormal subsegment γ| [0,t] . Notice that A ∪ {0} is a vector space, and its dimension is the corank of the abnormal path γ| [0,t] . Let λ be the unique normal lift of γ. Then for all α ∈ A ∪ {0} the family of curves

(4) s → γ α (s) = π • e (s-t) H (λ(t) + α), s ∈ [0, 1]
is a smooth family of normal geodesics, all coinciding with γ on the subinterval [0, t], and branching from it at time t. We know that each of the paths in this family is locally minimizing since they are normal. Moreover, if we take a compact subfamily of those, the time at which they are minimizing, starting from the branching point, can be chosen uniformly.

Theorem 7. Let γ α be a family of normal paths branching from γ at time t ∈ (0, 1), as in (4). Then for any compact subset A 0 ⊂ A ∪ {0} there exists ε > 0 such that γ a | [t-ε,t+ε] is the unique length-minimizing path between its extremities, up to reparametrizations, for all α ∈ A 0 .

The proof of Theorem 7 is a small adaptation of the proof for a single normal path, and it is an immediate consequence of the following more general result.

Proposition 8. Let x ∈ M , and let Λ ⊂ T *

x M be a compact set. For λ ∈ Λ, consider the normal paths γ λ (t) = π • e t H (λ). Then there exists ε > 0 such that, for all λ ∈ Λ, the restriction γ λ | [-ε,ε] is the unique length-minimizing path between its extremities, up to reparametrizations.

Proof. We want to apply the following obvious extension of [ABB19, Corollary 4.64].

Lemma 9. Let T > 0 and a ∈ C ∞ (M ). Let Ω 0 be an open subset of M such that, for all t ∈ [-T, T ], the map π • e t H • da| Ω 0 is a diffeomorphism from Ω 0 on its image Ω t . Let λ 0 ∈ L 0 ∩ π -1 (Ω 0 ) where L 0 = {d z a | z ∈ M } and define γ(t) = π • e t H (λ 0 ), for t ∈ [-T, T ].

Then γ is the unique length-minimizing path, up to reparametrization, among all horizontal paths γ : [-T, T ] → M with the same extremities and such that γ(t) ∈ Ω t for all t ∈ [-T, T ].

To do it, for λ ∈ Λ, we construct a family of functions a λ , continuous with respect to λ such that d x a λ = λ. Indeed, the theorem being a local result, we can suppose to be in a coordinates system (x 1 , . . . , x n ) on a neighborhood of x, and if λ = n i=1 λ i dx i , we define a λ (x 1 , . . . , x n ) = n i=1 λ i x i . Let Ω 0 be a relatively compact neighborhood of x which, for small T , contains

γ λ | [-T,T ] for all λ ∈ Λ. Consider the maps φ λ t = π • e t H • da λ | Ω 0 , for t ∈ [-T, T ] and λ ∈ T *
x M , noting that they are continuous in t and λ. For all λ ∈ Λ, we have φ λ 0 = Id| Ω 0 , so by semi-continuity of the rank, and by the fact that Ω 0 is compact, there exists a neighborhood of {0} × Λ where d x φ λ t in an isomorphism. By compactness of Λ, this neighborhood contains a set of the form [-t 0 , t 0 ] × Λ. By the inverse function theorem, and up to reducing Ω 0 , φ λ t is a diffeomorphism on its image for all (t, λ) ∈ [-t 0 , t 0 ] × Λ. Indeed, by compactness, the neighborhood of x given by the inverse function theorem can be uniformly chosen for (t, λ) ∈ [-t 0 , t 0 ] × Λ, by using a quantitative version of the latter, see [START_REF] Rifford | Sub-Riemannian geometry and optimal transport[END_REF]Theorem B.1.4].

Let K 1 ⊂ Ω 0 be a compact neighborhood of x. By continuity, there exists a neighborhood of {0} × Λ such that K 1 ⊂ Ω λ t = φ λ t (Ω 0 ) for all (t, λ) in this neighborhood. Since Λ is compact, we get t 1 ∈ (0, t 0 ] such that K 1 ⊂ Ω λ t for all t ∈ [-t 1 , t 1 ] and all λ ∈ Λ. Let then K 2 be a compact neighborhood of x included in the interior of K 1 and we find t 2 ∈ (0, t 1 ] such that γ λ (t) ∈ K 2 for all λ ∈ Λ and all t ∈ [-t 2 , t 2 ]. Finally, let us pose δ = d SR (K 2 , M \ K 1 ) > 0, and

ε = min t 2 , δ 4 2 max λ∈Λ H(λ) .
Let λ ∈ Λ and γ be a horizontal path defined for

[-ε, ε] such that γ(-ε) = γ λ (-ε) and γ(ε) = γ λ (ε), but whose image Γ is distinct from the image of γ λ . If Γ ⊂ K 1 , then γ(t)
∈ Ω λ t for all t, and we can thus apply Lemma 9 to conclude that (γ)

> (γ λ | [-ε,ε] ). Otherwise, there exists t * ∈ [-ε, ε] such that γ(t * ) / ∈ K 1 . Then, since γ(-ε) ∈ K 2 , we have: l(γ) ≥ δ > 2ε 2H(λ) = (γ λ | [-ε,ε] ).

An example of branching strictly normal geodesic

Let us stress that normal geodesics cannot branch in real-analytic sub-Riemannian structures, that is when the corresponding Hamiltonian function is real-analytic. In fact in this case, by the Cauchy-Kowalevski theorem, normal geodesics, which are projections of the solutions of the Hamiltonian equation, are real-analytic paths. By the principle of permanence, two distinct real-analytic paths cannot be equal on a segment. That is, the following well-known fact holds: Proposition 10. If H is real-analytic, normal geodesic cannot branch.

For building an example, we need to find a smooth, but non real-analytic structure, in which there is an abnormal geodesics that becomes strictly normal. A natural idea is to start from a structure admitting non-trivial abnormal geodesics (the simplest example being the flat Martinet structure) and "glue" it to a structure that do not admit non-trivial abnormal paths, like the Heisenberg structure. In fact, this works exactly as stated, and this is the idea that led us to the discovery of branching geodesics.

Let θ : R → [0, 1] be a smooth non-decreasing function such that θ(t) = 0 if t ≤ 0, θ(t) > 0 if t > 0, and θ(t) = 1 if t ≥ 1. Let A(x, y) = xθ(y) + x 2 θ(1 -y). Consider a rank 2 sub-Riemannian structure on R 3 defined by the following vector fields:

X = ∂ x Y = ∂ y + A(x, y)∂ z ,
so that we have a flat Martinet structure on the half-space y ≤ 0 and a Heisenberg one for y ≥ 1. The Lie bracket between those vector fields is:

[X, Y ] = B(x, y)∂ z , where B(x, y) = ∂ x A(x, y) = [θ(y) + 2xθ(1 -y)],
so that span{X, Y, [X, Y ]} = R 3 , except on the so-called Martinet surface

Σ = {(x, y, z) ∈ R 3 | B(x, y) = 0}.
At all points in Σ, we have [X, [X, Y ]] = 2θ(1 -y)∂ z = 0, therefore Σ is smooth and the distribution is bracket-generating. It is well-known that abnormal paths for this distribution are exactly the horizontal ones contained in Σ, see for example [START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF][START_REF] Rifford | Sub-Riemannian geometry and optimal transport[END_REF]. To characterize normal geodesics, the Hamiltonian function is

H = 1 2 h 2 X 1 + h 2 X 2 = 1 2 p 2 x + (p y + A(x, y)p z ) 2 .
The Hamiltonian vector field is thus, in coordinates (x, y, z, p x , p y , p z ):

(5)

H =              ∂H ∂px ∂H ∂py ∂H ∂pz -∂H ∂x -∂H ∂y -∂H ∂z              =             p x p y + A(x, y)p z (p y + A(x, y)p z )A(x, y) -p z (p y + A(x, y)p z )B(x, y) -p z (p y + A(x, y)p z )∂ y A(x, y) 0             .
In particular, the path in the cotangent bundle (0, t, 0, 0, 1, 0), for t ∈ R is an integral curve of H, and therefore the lift of the normal geodesic γ(t) = (0, t, 0). For t < 0, its projection γ is contained in the Martinet surface Σ, and therefore this part of the curve is abnormal. Indeed, for every α = 0, the path (0, t, 0, 0, 0, α) is an abnormal lift of this geodesic. As soon as t > 0, however, γ is not contained in Σ and therefore any such a segment is strictly normal. It is quite natural for this to happen as the Heisenberg structure has no abnormal geodesics. So if we consider the path γ(t) = (0, t, 0) for t ∈ [-T, T ] for some T > 0, it has a maximal abnormal subsegment [-T, 0] and therefore by Corollary 6, it branches at time t = 0. What happens is that, starting at t = -T and until t = 0, the differential of the end-point map has a 1-dimensional cokernel for the corresponding control, which is the family of covectors (0, 0, α) for α ∈ R, and thus the space of initial covectors of normal lifts for this path is the 1-dimensional affine space {(0, 1, α) | α ∈ R}. Once time t = 0 is attained the abnormal geodesic can still be prolonged (as the trace of the distribution in Σ) but it loses its normal status, becoming strictly abnormal. Meanwhile, the 1-dimensional family of normal lifts can be prolonged yielding a family of distinct geodesics, which are all strictly normal. In Figure 1 we computed, using the Euler method, some of those geodesics γ α , for different values of α.

Finally, we observe that the collection of those normal geodesics do describe an embedded surface of R 3 , at least locally around the y-axis (which is the normal geodesic with initial covector (0, 1, 0)) as shown in this result: Proposition 11. The map Φ : R 2 → R 3 defined by Φ(t, α) = π • e (t+T ) H (λ α ), where λ α is the initial covector (0, 1, α) at point (0, -T, 0), is an embedding on a neighborhood of any point (t, 0) with t > 0.

Proof. Since Φ(t, 0) = (0, t, 0), we have ∂Φ ∂t (t, 0) = ∂ y . So we just have to show that ∂Φ ∂α (t, 0) is independent from ∂ y , which we will do by pointing out that its x component is non zero. Let t → (x α (t), y α (t), z α (t), p α x (t), p α y (t), p α z (t))

Figure 1. Numerical plot of the branching geodesics γ α , projected on the xy plane. Notice that the abnormal path lies in the Martinet which must bend in order to avoid the Heisenberg region.

be the solution of Hamilton's equation ( 5), with initial condition (0, -T, 0, 0, 1, α) at time t = 0, in such a way that Φ(t, α) = (x α (t), y α (t), z α (t)). We observe that p α z (t) = α for all times, and thus:

x α (t) = t 0 p α x (s)ds = t 0 s 0 ṗα x (τ )dτ ds = -α t 0 s 0 f (α, τ )dτ ds.
Here f (α, τ ) is a smooth function whose expression can be obtained from (5). We notice only that f (0, τ ) = θ(τ ). So

∂x α (t) ∂α α=0 = - t 0 s 0 θ(τ )dτ ds,
which is non-zero by the properties of θ.

Remark 12. As anticipated in Section 1.1, any smooth A such that A(0, y) = 0, and such that B = ∂ x A satisfies B(x, y) = x for y < 0, and B(0, y) > 0 when y > 0, yields a one-parameter family of branching geodesics, verifying Proposition 11.

Normal geodesics with multiple branching

In the example from the previous section, the corank function of a path γ α starting at t = -T is equal to one for t ≤ 0 and zero for t > 0. From this, we can easily construct any kind of corank function, and therefore any kind of normal branching. To do that, we consider a suitable product of our example.

The product M 1 × M 2 of sub-Riemannian manifolds M 1 and M 2 has a product sub-Riemannian structure simply defined as the direct sum of the distributions and the metrics in M 1 and M 2 . It is then easy to see that a path γ is a geodesic in M if and only if γ 1 and γ 2 are geodesics in M 1 and M 2 , where γ i is the projection of γ on M i . Moreover, if λ 1 and λ 2 are normal Lagrange multipliers for γ 1 and γ 2 respectively, (λ 1 , λ 2 ) is a normal Lagrange multiplier for γ, and the converse is true. The same is true for abnormal multipliers, so the corank of a path γ is the sum of the coranks of the γ i . From these observations, we get the following result, which gives a complete description of what kind of branching one can expect: Theorem 13. Let f : [0, 1] → N be a nonincreasing left-continuous function. Then there exists a sub-Riemannian manifold M and a normal geodesic γ on M such that its corank function coincides with f . In particular γ branches at the jumps of f . Proof. Constant paths have corank equal to the corank of the distribution, therefore the corank of M must be equal to f (0). So we define M as the sub-Riemannian product M = N f (0) , where N is the sub-Riemannian structure in R 3 defined in Section 4. Denote then by a = f (0) -f (0 + ) (f (0 + ) being the right limit of f at 0), b = f (1) and t 1 , . . . , t k the times of the discontinuities of f in (0, 1), each one repeated multiple times according to the amplitude of the discontinuity. Denote as previously by γ the path in N defined by γ(t) = (0, t, 0) for t ∈ R and take γ 0 any strictly normal geodesic with no non-trivial abnormal subsegments in N (for example γ 0 = γ| [0,1] ).

If we define the path γ : [0, 1] → M by: γ(t) = (γ 0 (t), . . . , γ 0 (t) a times

, γ(t -t 1 ), γ(t -t 2 ), . . , γ(t -t k ), 0 R 3 , . . . , 0 R 3 b times

), then the corank function of γ is exactly f . As a remark, we also have a full description of the normal geodesics branching from γ, given, for α 1 , . . . , α k ∈ R, by: γ(α 1 ,...,α k ) (t) = (γ 0 (t), . . . , γ 0 (t), γ α 1 (t -t 1 ), γ α 2 (t -t 2 ), . . . , γ α k (t -t k ), 0 R 3 , . . . , 0 R 3 ), where γ α are the already defined branching geodesics in N , from Section 4.

Acknowledgments. This work was supported by the Grants ANR-15-CE40-0018 and ANR-15-IDEX-02. While this note was in preparation, N. Juillet shared with us an example similar to the one of Section 4, in a rank-varying structure on R 4 . We thank the anonymous referee for carefully reading the paper, and suggesting the simple description of branching in terms of magnetic field presented in Section 1.1.