Pierre-Louis Giscard

Exact solutions for the time-evolution of quantum spin systems under arbitrary waveforms using algebraic graph theory

Keywords: exact solutions, time-ordered exponentials, time evolution, quantum spin systems, path sum, algebraic graph theory

A general approach is presented that offers exact analytical solutions for the time-evolution of quantum spin systems during parametric waveforms of arbitrary functions of time. The proposed method utilises the path-sum method that relies on the algebraic and combinatorial properties of walks on graphs. A full mathematical treatment of the proposed formalism is presented, accompanied by an implementation in Matlab. Using computation of the spin dynamics of monopartite, bipartite, and tripartite quantum spin systems under chirped pulses as exemplar parametric waveforms, it is demonstrated that the proposed method consistently outperforms conventional numerical methods, including ODE integrators and piecewise-constant propagator approximations.

Introduction

Understanding the dynamics and control of quantum spin systems, indispensable in spectroscopy, sensing, quantum computing and information processing, is among the most challenging areas of research in current science. The diverse fields of applications include high-resolution magnetic resonance spectroscopy and imaging [START_REF] Khaneja | Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms[END_REF], terahertz technologies [START_REF] Rasanen | Optimal control of quantum rings by terahertz laser pulses[END_REF][START_REF] Coudert | Optimal control of the orientation and alignment of an asymmetric-top molecule with terahertz and laser pulses[END_REF], and control of trapped ions [START_REF] Zhao | Coherent and optimal control of adiabatic motion of ions in a trap[END_REF], cold atoms [START_REF] Saywell | Optimal control of mirror pulses for cold-atom interferometry[END_REF] and NV-centers in diamond [START_REF] Chou | Optimal control of fast and high-fidelity quantum gates with electron and nuclear spins of a nitrogen-vacancy center in diamond[END_REF][START_REF] Tian | Optimal quantum optical control of spin in diamond[END_REF]. In many of these applications, sophisticated manipulations of quantum systems are achieved using pulses in the form of radio-frequency, microwave, or laser. A good understanding of the evolution of quantum spin systems during these events is crucial for enabling new methodologies.

A large class of these pulses are parametric, i.e. the waveform can be represented as a function dependent on certain time-dependent and time independent parameters. Numerical solutions of the spin dynamics during such pulses can be presented via solutions of ODEs e.g. using adaptive Runge-Kutta method [START_REF] Foroozandeh | Spin dynamics during chirped pulses: applications to homonuclear decoupling and broadband excitation[END_REF], or methods relying on the approximation of matrix exponentials and Fokker-Planck formalism [START_REF] Kuprov | Fokker-planck formalism in magnetic resonance simulations[END_REF][START_REF] Allami | Quantum mechanical mri simulations: Solving the matrix dimension problem[END_REF], but exact solutions of quantum systems under arbitrary parametric pulses are less explored. Although for frequency-swept pulses, an important class of parametric pulses, some approximated [START_REF] Dumez | Spatial encoding and spatial selection methods in highresolution NMR spectroscopy[END_REF] and exact [START_REF] Hioe | Solution of bloch equations involving amplitude and frequency modulations[END_REF][START_REF] Zhang | Full analytical solution of the bloch equation when using a hyperbolic-secant driving function[END_REF] analytical solutions for the time evolution of single spin-1 {2 particles have been presented in the literature. Other Email addresses: giscard@univ-littoral.fr (Pierre-Louis Giscard), mohammadali.foroozandeh@chem.ox.ac.uk (Mohammadali Foroozandeh) approaches that can offer similar solutions include integrable multi-state Landau-Zener Models [START_REF] Sinitsyn | The quest for solvable multistate landauzener models[END_REF][START_REF] Sinitsyn | Integrable time-dependent quantum hamiltonians[END_REF] In this article we present very general exact analytical solutions for spin dynamics during parametric waveforms of arbitrary functions of time using the path-sum approach [START_REF] Giscard | Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR[END_REF][START_REF] Giscard | An exact formulation of the time-ordered exponential using path-sums[END_REF]. The approach relies on the algebraic and combinatorial properties of walks on graphs and here yields the complete description of the time-dependent propagator matrix in terms of the waveform and spin system parameters. The proposed approach is general and applicable to any user-defined pulses, arbitrarily constructed of parametric components e.g. chirps, hyperbolic secants, polynomials, Fourier series, etc. The method can be utilised to describe the time-evolution of quantum systems under the effect of parametric pulses. This paper is structured as follows: in section 2 the underlying theory for the equation of motion of quantum spin systems along with a general background on the path-sum approach is presented. In section 3 we present in detail the analytical solution for the propagation of monopartite systems driven by arbitrary time-dependent pulses in the path-sum formalism. In section 4 the method is extended to larger systems of spin-1 {2 particles and explicitly worked out in the case of bipartite and tripartite systems, demonstrating the applicability of the path-sum formalism to diverse Hamiltonian structures. Finally, in section 5 we present numerical results and assess the computational performances of Matlab codes implementing the path-sum formalism, using 6-parameter chirped pulses as exemplar waveforms, in the context of nuclear magnetic resonance (NMR) and the presence of scalar interactions.

Theory

2.1. Quantum equation of motion for a single spin- 1 2 Generators of rotation in SOp3q can be expressed as orthogonal skew-symmetric matrices

L x " ¨0 0 0 0 0 ´1 0 1 0 ', L y " ¨0 0 1 0 0 0 ´1 0 0 ', L z " ¨0 ´1 0 1 0 0 0 0 0 ', (1)
and the Hamiltonian of the system under a pulse can be written as H ptq "

2β x ptqL x `2β y ptqL y `ΩL z , (2)
where Ω is the spin resonance offset, as commonly used in nuclear magnetic resonance (NMR), and β x and β y are real and imaginary components of the pulse

βptq " β x ptq `iβ y ptq " 1 2 ω 1 ptq exp `iφptq ˘. (3)
The state of the system, ρ can be expressed as

ρptq " g 1 ptqL x `g2 ptqL y `g3 ptqL z , (4)
where gptq " rg 1 ptq, g 2 ptq, g 3 ptqs P R 3 is a time-dependent unit vector representing the position of the spin on a Bloch sphere. Considering 9 gptq " ´iH ptq gptq,

we can write

¨9 g 1 9 g 2 9 g 3 '" ¨0 ´Ω 2β y Ω 0 ´2β x ´2β y 2β x 0 '¨g 1 g 2 g 3 '. (6)
Here time-dependencies of g 1 , g 2 , g 3 , β x , and β y are removed for simplicity. For a single spin-1 2 the basis set can also be written using shift operators L "

! L `, ? 2L z , L ´) , (7)
where L `" L x `iL y and L ´" L x ´iL y . Using this basis set, the Hamiltonian under an arbitrary waveform can be written as

H ptq " βptqL ``βptqL ´`ΩL z , (8)
where the bar ¯indicates complex conjugation. Equation (5) using this basis set can be written as

¨9 g 1 9 g 2 9 g 3 '" ´i ¨´Ω ? 2β 0 ? 2 β 0 ´?2β 0 ´?2 β Ω '¨g 1 g 2 g 3 '. (9)
Switching between two representations of the Hamiltonian in Equations ( 6) and ( 9) can be achieved with a transformation matrix U as

H S Z " U ¨H C ¨U: , (10)
where

U " 1 ? 2
¨1 i 0 0 0 ? 2 1 ´i 0 ', (11)
and H S Z and H C are the Hamiltonians in the shift and Cartesian basis respectively.

Path-sum approach

The system eq. ( 9) is solved exactly with the method of pathsum, which we here present only briefly. We refer the reader to [START_REF] Giscard | Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR[END_REF][START_REF] Giscard | An exact formulation of the time-ordered exponential using path-sums[END_REF] for further details. Conceptually, the method relies on three mathematical pillars:

i) The evolution operator solution of the quantum equation of motion is the resolvent of the Hamiltonian H if the usual multiplication between entries of H is replaced with another product, called Volterra composition [START_REF] Volterra | Lec ¸ons sur la composition et les fonctions permutables, Collection de monographies sur la théorie des fonctions[END_REF].

ii) This resolvent is formally given by a sum over all walks on a graph G H whose adjacency matrix is the Hamiltonian operator and which thus encodes the discrete structure of the quantum state space.

iii) Sums of walks are given exactly by a continued fraction of finite depth and breadth involving a finite number of progressively simpler terms and stemming from an algebraic structure of sets of walks. Remarkably, each term of this continued fraction corresponds to a simple cycle or a simple path on the graph, i.e. to walks that do not visit any vertex more than once [START_REF] Giscard | Walk-sums, continued fractions and unique factorisation on digraphs[END_REF].

Combining these three observations, one can evaluate any time-ordered evolution operator exactly from a finite number of Volterra compositions and inverses between entries of the Hamiltonian as they appear successively along simple cycles and simple paths on the quantum state space. This conceptual framework has already been established in its full generality, so that when implementing path-sum concretely only writing the continued fraction needs to be done. This is because it changes according to the situation since the structure of the Hamiltonian dictates what the quantum state space graph G H looks like and thus which simple cycles/paths occur on it.

We begin by recasting eq. ( 9) in matrix form as d dt U " ´iH ptqU ptq with H the Hamiltonian matrix of eq. ( 9) and U the corresponding evolution operator, which is the time-ordered exponential of ´iH ,

U ptq " T exp " ´i ż t 0 H pt 1 q dt 1  , (12)
As stated earlier, we can put this in resolvent form upon using the Volterra composition, a special case of the ‹-product [START_REF] Giscard | On the solutions of linear Volterra equations of the second kind with sum kernels[END_REF].

-Ω

2 β(t) 2 β(t) -2 β(t) -2 β(t) Ω 1 2 3
Figure 1: Graph G H showing the structure of the quantum state space as imposed by the Hamiltonian H . Edge weights are indicated next to each edge while the vertices are labelled by integers 1, 2 and 3. The vertices in fact correspond to the canonical basis states of the quantum state space, 1 " |1, 0, 0y, 2 " |0, 1, 0y and 3 " |0, 0, 1y. Alternative bases are of-course possible, as path-sum is valid in all bases. The contribution of a walk w on G H to the overall quantum dynamics as described by U is given the ‹-products of the weights of the edges it traverses. This allows for an exact evaluation of U using a path-sum which represents the sum of all weighted walks on the above graph.

remark 1. Fundamentally, the ‹-product is defined differently and on a much wider class of distributions of two variables, which include piecewise-smooth functions but also the Dirac delta and all its distributional derivatives. From there, it is also defined on matrices of distributions [START_REF] Giscard | On the solutions of linear Volterra equations of the second kind with sum kernels[END_REF]. These technicalities are necessary to rigorously prove, among other things, that the Dirac delta distribution is the identity element with respect to the ‹-product and to perform ‹-inversions. Here we admit such results and proceed with the simplest definition of the ‹-product as Volterra did before the inception of distributions [START_REF] Volterra | Lec ¸ons sur la composition et les fonctions permutables, Collection de monographies sur la théorie des fonctions[END_REF].

To define the resolvent form of U ptq, let f pt 1 , tq and gpt 1 , tq be two smooth functions of two variables. We define f ‹ g as

`f ‹ g ˘pt 1 , tq " ż t 1 t f pt 1 , τqgpτ, tqdτ, (13)
For functions of less than two time variables, the variable must be treated as the left one, e.g. if hpt 1 q depends only on one variable, then ph ‹ gqpt 1 , tq " hpt 1 q ş t 1 t gpτ, tqdτ, p f ‹ hqpt 1 , tq " ş t 1 t f pt 1 , τqhpτqdτ. One can then show that [START_REF] Giscard | An exact formulation of the time-ordered exponential using path-sums[END_REF] U ptq "

ż t 0 ´I‹ ´p´iqH ¯‹´1 pτ, 0q dτ,
where I ‹ :" 1 ‹ ˆI " δpt 1 ´tqI is the identity matrix times the Dirac delta distribution 1 ‹ " δpt 1 ´tq and p.q ‹´1 designates a matrix inverse where all ordinary multiplications are taken instead to be Volterra compositions. To alleviate the notation and because it is related to Green's functions, the above resolvent p1 ‹ I ´H q ‹´1 pt 1 , tq is denoted G pt 1 , tq. Then U ptq " ş t 0 G pτ, 0qdτ. Now path-sum guarantees not only that G pt 1 , tq is well defined, but also that any entry of G is given from ‹-products and ‹-inverses of entries ´iH . This is illustrated in details in Section 3.

Analytical solution of the quantum equation of motion

Solutions for general parametric waveforms

We consider the Hamiltonian of eq. ( 8), with the corresponding graph G H illustrated in Figure 1. The path-sum formula- ¯‹´1 .

Here all ´i factors originates from the ´i in ´iH . The ‹inverses in the above expression, e.g. in p1 ‹ ´p´iqΩq ‹´1 , stem from sums over all possible repetitions of a simple cycle c. Indeed, the observation 1 `c `c.c `c.c.c `¨¨¨" ř n c n " p1 ´cq ´1 is true at the formal level for cycles and walks [START_REF] Giscard | Walk-sums, continued fractions and unique factorisation on digraphs[END_REF], which implies its validity for actual weighted walks [START_REF] Giscard | An exact formulation of the time-ordered exponential using path-sums[END_REF].

We can now evaluate the above path-sum explicitly. First, since Ω depends on strictly less than two time variables we have [START_REF] Giscard | Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR[END_REF] `1‹ ´p´iqΩ ˘‹´1 pt 1 , tq " δpt 1 ´tq ´iΩe ´iΩpt 1 ´tq .

Second, since β depends on a single time, for any function f pt 1 , tq we have

pβ ‹ f qpt 1 , tq " βpt 1 q ż t 1 t f pτ, tqdτ.
With these observations, we obtain

G 22 pt 1 , tq " ˆ1‹ ´p´2qβpt 1 q ż t 1 t βpτqdτ ´p´2q βpt 1 q ż t 1 t βpτqdτ ´p´2q βpt 1 q ż t 1 t `eiΩpt 1 ´τq ´1˘β pτqdτ ´p´2qβpt 1 q ż t 1 t `e´iΩpt 1 ´τq ´1˘β pτqdτ ˙‹´1 " ˆ1‹ `2 B Bt 1 ˇˇż t 1 t e ´iΩτ βpτqdτ ˇˇ2 ˙‹´1 , " ˆ1‹ `2 B Bt 1 ˇˇB t 1 ,t pΩ{2πq ˇˇ2 ˙‹´1 . (14)
In these expressions, B t 1 ,t pΩ{2πq :" ş t 1 t e ´iΩτ βpτqdτ is the Fourier transform with respect to τ of βpτq `Θpτ´t 1 q´Θpτ´tq ȇvaluated in Ω{2π, Θp.q being the Heaviside function with the convention that Θp0q " 1. remark 2. Appearance of Fourier transform in the path-sum solution in eq. ( 14) is not a feature of the path-sum method itself, but rather of the peculiar form of the Hamiltonian H , in particular that Ω is time-independent.

The end result is remarkably simple and holds for all pulse shapes β, provided they are smooth in the mathematical sense. Numerical evaluation of the exact path-sum solution eq. ( 14) is straightforward using a discretized version of the ‹-product. Analytically speaking, the path-sum solution is best evaluated from its Neumann expansion, a non-perturbative series representation that is super-exponentially (and thus unconditionally) convergent and lends itself to analytical calculations. In order to facilitate the notation, we introduce bpt 1 , tq :"

B Bt 1 ˇˇB t 1 ,t pΩ{2πq ˇˇ2, " 2 Re ˜e´iΩt 1 βpt 1 q ż t 1 t e iΩτ βpτqdτ ¸. (15)
Then, since

U 22 ptq " 1 `ż t 0 G 22 pτ, 0qdτ, (16)
the Neumann expansion of the path-sum solution is

U 22 ptq " 1 `8 ÿ n"1 p´2q n ż t 0 b ‹n pτ, 0qdτ, " 1 ´2 ˇˇB t,0 pΩ{2πq ˇˇ2 (17) `4 ż t 0 ˇˇB t,τ pΩ{2πq ˇˇ2 B Bτ ˇˇB τ,0 pΩ{2πq ˇˇ2 dτ ´¨¨¨,
here displaying only the first two orders. Higher order terms are analytically accessible and can reach any desired accuracy. To illustrate this, consider the m th Neumann order approximation

U pmq 22 ptq :" 1 `m ÿ n"1 p´2q n ż t 0 b ‹n pτ, 0qdτ. (18)
In this notation U p8q 22 ptq is the exact solution, which we emphasize that it can be evaluated numerically directly without resorting to Neumann series. Indeed, Volterra compositions map to ordinary matrix products when time is discretized [START_REF] Birkandan | Computations of general Heun functions from their integral series representations[END_REF] so that the ‹-inverse of eq. ( 14) becomes the ordinary inverse of a well-conditioned triangular matrix, directly yielding the numerical evaluation of U Using path-sum we can evaluate the other entries of G similarly.

G 11 ptq " ´1‹ ´Loop 1Ð1 hk kik kj iΩ `2 β ‹ `1‹ `Cycle 2Ð3Ð2 hkkkkkkkkkkkkkkkikkkkkkkkkkkkkkkj 2β ‹ p1 ‹ `iΩ lo omo on Loop 3Ð3 q ‹´1 ‹ β ˘‹´1 ‹ β looooooooooooooooooooooooooooomooooooooooooooooooooooooooooon Cycle 1Ð2Ð1 ¯‹´1 , " ´1‹ ´iΩ`2β‹ ´1‹ `2e ´iΩt 1 βpt 1 q Ě B t 1 ,t pΩ{2πq ¯‹´1 ‹ β¯‹ ´1,
and G 33 ptq " Ğ G 11 ptq. Off-diagonal terms are given by sums over simple paths on the graph of Figure 1. Here we have for example

G
where U ptq is the solution of the differential equation

d dt U ptq " ´iH ptqU ptq, U p0q " 1. (20)
In general, regardless of whether the Hamiltonian commutes with itself or not, the solution is eq. ( 12).

Solution for monopartite systems

The Hamiltonian of a monopartite quantum spin-1 {2 system with a resonance offset Ω can be written as

H ptq " ˆΩ 2 βptq βptq ´Ω 2 ˙. (21)
Since the Hamiltonian is 2 ˆ2 with no zero entry, the corresponding graph G H is the complete graph on 2 vertices and the path-sum formulation of the corresponding evolution operator is [START_REF] Giscard | Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR[END_REF]

G 11 pt 1 , tq " ´1‹ `Loop 1Ð1 hk kik kj iΩ{2 `Edge 1Ð2 hk kik kj β ‹ p1 ‹ ´Loop 2Ð2 hk kik kj iΩ{2 q ‹´1 ‹ Edge 2Ð1 hk kik kj β looooooooooooooooooooooomooooooooooooooooooooooon Cycle 1Ð2Ð1 ¯‹´1 , " ´1‹ `iΩ{2 `βpt 1 qe iΩt 1 {2 B t 1 ,

System of M interacting spin-1 2 s

In an M-spin system, the general form of spin operator for the i th spin can be expressed as (22) with a Pauli matrix in position i. The total Hamiltonian can be written as a sum of the Hamiltonian for spin resonance offset (H Ω ), pulse (H P ), and some scalar interaction, known as Jcoupling in NMR (H J )

I piq α " 1 b 1 b ¨¨¨b σ α b ¨¨¨1 b 1 α P tx, y, zu
H Ω " M ÿ i"1 Ω i I piq z , H P " f ptq M ÿ i"1 I piq x `gptq M ÿ i"1 I piq y , H J " M ÿ i, j jąi

Solution for bipartite systems

The Hamiltonian for a bipartite quantum spin-

where

h 11 " 1 2 pπJ `Ω1 `Ω2 q h 22 " 1 2 p´πJ `Ω1 ´Ω2 q h 33 " 1 2 p´πJ ´Ω1 `Ω2 q h 44 " 1 2 pπJ ´Ω1 ´Ω2 q
For the sake of simplicity regarding the equations that give the evolution operator, and without loss of generality, let us change the energy gauge to

H 1 ptq " H ptq ´1 2 πJ I ,
with I the 4 ˆ4 identity matrix. Of course, this Hamiltonian leads to the very same dynamics as of the original H of eq. ( 23). Now we rely on path-sum's scale invariance to express the solution of this spins system. Mathematically speaking, this 'scale invariance' refers to path-sum's continued validity for all partitions of the Hamiltonian matrix into (possibly noncontiguous) blocks. In other terms, there is a path-sum formulation of the evolution operator U ptq in terms of any chosen partition of the Hamiltonian into blocks of any size, including non-square ones. Physically, this means that path-sum is able to formulate the dynamics of a quantum system in terms of the isolated dynamics of any chosen collection of its subsystems. We may also use partitions exhibiting mappings from one Hamiltonian to another, so as to demonstrate the similarity in the time evolutions they effect.

To illustrate the direct relation between monopartite and bipartite dynamics, we partition the Hamiltonian H 1 ptq as follows

H 1 ptq " ¨h1 11 βptqu 0 βptqu T H 1 II βptqu T 0 βptqu ´h1 11 ', (24)
where we defined u :" p1, 1q and

H 1 II :" ˆh1 22 πJ πJ h 1 33 ˙,
is the 2 ˆ2 submatrix of H ptq formed by entries H 1 ptq i, j with i, j " 2, 3. In these expressions, h 1 ii " h ii ´p1{2qπJ. We chose this partition so has to exhibit an exact mapping from the bipartite system to an effective monopartite case with non-Abelian energy H II for the degrees of freedom spanned by states p0, 1, 0, 0q and p0, 0, 1, 0q. This mapping is best seen on the graph G H 1 , shown in Figure 2, which encodes this partition of H 1 and is structurally similar to that of one the spin case shown in Figure 1. The partition in itself plays no special role in the analytical solution however, and the problem could be solved using the ordinary entries of H 1 directly. Now let U ptq be the evolution operator of the bipartite system-i.e. the solution of [START_REF] Giscard | On the solutions of linear Volterra equations of the second kind with sum kernels[END_REF] with Hamiltonian H 1 -and let G " 9 U . Using the same partition for U as for H 1 , let U II be the 2 ˆ2 submatrix of U formed by the entries U ptq i, j with h' 11

β u β u T β u T β u -h' 11 H ' II 1 II 4 Figure 2: Graph G H 1
showing the structure of the quantum state space as imposed by the bipartite Hamiltonian H 1 when partitioned as per eq. ( 24).

Because the structure of G H 1 and of the monopartite Hamiltonian graph G H of Figure 1 differ only in the presence of a central loop on vertex 2, the pathsum formulation of the corresponding evolution operators will differ only in a single term representing this loop.

i, j " 2, 3. Then U II ptq " ş t 0 G II pτqdτ with G II the Green's function given by the path-sum expression

G II pt 1 q " ˆI‹ `iH 1 II `2P B Bt 1 ˇˇB t 1 ,t `pΩ 1 `Ω2 q{4π ˘ˇ2 ˙‹´1 ,
where I ‹ " 1 ‹ ˆI " δpt 1 ´tqI is the two by two identity matrix times 1 ‹ " δpt 1 ´tq and P " u T .u is the 2 ˆ2 matrix full of 1. As predicted above, G II is similar to [START_REF] Sinitsyn | The quest for solvable multistate landauzener models[END_REF] for G 22 in the monopartite case. Note that the additional iH 1 II ptq term would also have been present in the one spin case had there been a non-zero H 22 diagonal term in the Hamiltonian. The pathsum for G II involves a matrix ‹-inverse, meaning that G II enters a non-Abelian linear Volterra integral equation of the second type, with unconditionally convergent analytical representation given by the matrix-valued Neumann series of ‹-powers of its kernel.

Furthermore, one can use path-sum's scale invariance to get scalar-valued expressions for any entry of G II rather than manipulating matrices. For example, entry pG II q 1,1 " G 2,2 is given by the path-sum The remaining entries are succinctly obtained from path-sum in matrix form as

pG II q 1,1 " ˆ1‹ `ih 1 22 `2bpt 1 , tq ´p2bpt 1 , tq `iπJq ‹ `1‹ `2bpt 1 , tq `ih 1 33 ˘‹´1 ‹ p2bpt 1 , tq `iπJq ˙‹´1 , where bpt 1 , tq " B Bt 1 ˇˇB t 1 ,t `pΩ 1 `Ω2 q{4π ˘ˇ2 . Vector G 1,II " pG 1,2 , G 1,3 q then follows as G 1,II " pG 1,2 , G 1,3 q, " p1 ‹ ´ih 1 11 q ‹´1 ‹ p´i βuq ‹ G II ,
ˆG1 G 1,4 G 4,1 G 4 ˙" ˜I‹ ´ˆ´ih 1 11 0 0 ih 1 11 βpt 1 q βpt 1 q βpt 1 q βpt 1 q ˙.

Solution for tripartite systems

The Hamiltonian for a tripartite quantum spin-1 {2 system with offsets Ω 1 , Ω 2 , and Ω 3 and coupling constants J 12 , J 13 , and J 23 can be expressed as

H ptq " (25
β 0 0 0 β 0 β β h 88 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ '
, where h 11 " 1 2 pπ pJ 12 `J13 `J23 q `Ω1 `Ω2 `Ω3 q , h 22 " 1 2 pπJ 12 ´π pJ 13 `J23 q `Ω1 `Ω2 ´Ω3 q , h 33 " 1 2 p´π pJ 12 ´J13 `J23 q `Ω1 ´Ω2 `Ω3 q , h 44 " 1 2 p´π pJ 12 `J13 ´J23 q `Ω1 ´Ω2 ´Ω3 q , h 55 " 1 2 p´π pJ 12 `J13 ´J23 q ´Ω1 `Ω2 `Ω3 q , h 66 " 1 2 p´π pJ 12 ´J13 `J23 q ´Ω1 `Ω2 ´Ω3 q , h 77 " 1 2 pπJ 12 ´π pJ 13 `J23 q ´Ω1 ´Ω2 `Ω3 q , h 88 " 1 2 pπ pJ 12 `J13 `J23 q ´Ω1 ´Ω2 ´Ω3 q .

(26

)
This system presents no further difficulty than the monopartite and bipartite cases. We may exploit different partitions of the quantum state space to exhibit mappings either from the tripartite to the bipartite or monopartite cases, or to an altogether

h 11 h 88 β(t) u T β(t) u β(t) u T β(t) u β(t) M β(t) M 1 8
H II H III

II III

Figure 3: Graph G H showing the structure of the quantum state space as imposed by the tripartite Hamiltonian H when partitioned as per Eq. ( 27). Different partitions of the Hamiltonian of eq. ( 25) would lead to different graphs, all of which generate valid path-sum expressions of the very same solution.

different situation such as a mathematically advantageous structure, all thanks to path-sum scale invariance. Here we choose to map the Hamiltonian to a path-graph, whose path-sums yield continued fractions with a single branch.

Let |iy designate the canonical basis states, e.g. |1y " p1, 0, 0, 0, 0, 0, 0, 0q T .

Let V I " spanp|1yq, V II " spanp|2y, |3y, |5yq, V III " spanp|4y, |6y, |7yq and V IV " spanp|8yq. At the scale formed by these vector spaces, the Hamiltonian takes the form

H ptq " ¨h11 βptqu 0 0 βptqu T H II βptqM 0 0 βptqM H III βptqu T 0 0 βptqu h 88 ‹ ‹ ' , (27)
where u " p1, 1, 1q and

M " ¨1 1 0 1 0 1 0 1 1 ',

'.

The corresponding graph is a path-graph on 4 vertices, illustrated on Figure 3.

To succinctly present the path-sum expression for the evolution operator let G " 9 U , define I ‹ :" I ˆ1‹ the 3-by-3 identity matrix times the Dirac Delta distribution 1 ‹ " δpt 1 ´tq and let

Γ 8 " 1 1 ‹ ´p´iqh 88
" δpt 1 ´tq ´ih 88 e ´ih 88 pt 1 ´tq ,

Γ III " 1 I ‹ ´p´iqH III ´p´iq 2 βu T ‹ Γ 8 ‹ βu , Γ II " 1 I ‹ ´p´iqH II ´p´iq 2 βM ‹ Γ III ‹ βM , γ 1 " 1 1 ‹ ´p´iqh 11 " δpt 1 ´tq ´ih 11 e ´ih 11 pt 1 ´tq , γ II " 1 I ‹ ´p´iqH II ´p´iq 2 βu T ‹ γ 1 ‹ βu , γ III " 1 I ‹ ´p´iqH III ´p´iq 2 βM ‹ γ II ‹ βM ,
here all inverses are to be understood as ‹-inverses, i.e. inverses with respect to the ‹-product, the above presentation being chosen so as to reveal clearly the continued fraction nature the pathsum expressions. Now we have access to all entries and block of entries of the Green's function G ,

G " ¨G1,1 G 1,II G 1,III G 1,8 G II,1 G II,II G II,III G II,8 G III,1 G III,II G III,III G III,8 G 1,8 G II,8 G III,8 G 8,8 ‹ ‹ ' ,
with the path-sum expressions

G 1,1 " `1‹ ´p´iqh 11 ´p´iq 2 βu ‹ Γ II ‹ βu T ˘‹´1 , G II,II " ´I‹ ´p´iqH II ´p´iq 2 βu T ‹ γ 1 ‹ βu ´p´iq 2 βM ‹ Γ III ‹ βM ¯‹´1 , G III,III " ´I‹ ´p´iqH III ´p´iq 2 βu T ‹ Γ 8 ‹ βu ´p´iq 2 βM ‹ γ II ‹ βM ¯‹´1 , G 8,8 " `1‹ ´p´iqh 88 ´p´iq 2 βu ‹ γ III ‹ βu T ˘‹´1 ,
while the off-diagonal blocs are

G II,1 " Γ II ‹ βu T ‹ G 1,1 , G III,8 " γ III ‹ βu T ‹ G 8,8 , G III,1 " Γ III ‹ βM ‹ G II,1 , G II,8 " γ II ‹ βM ‹ G III,8 , G 8,1 " Γ 8 ‹ βu ‹ G III,1 , G 1,8 " γ 1 ‹ βu ‹ G II,8 , and
G 1,II " γ 1 ‹ βu ‹ G II,II , G 8,III " Γ 8 ‹ βu ‹ G III,III , G III,II " Γ III ‹ βM ‹ G II,II , G II,III " γ II ‹ βM ‹ G III,III , G 8,II " Γ 8 ‹ βu ‹ G III,II , G 1,III " γ 1 ‹ βu ‹ G II,III .
This provides the exact, fully analytical formulation of the evolution of the tripartite systems under a time-dependent driving.

These expressions can all be expanded analytically into closed form approximations via truncated unconditionally convergent Neumann series. They can also be directly evaluated numerically when time is discretized using strategies presented in the following Section.

Numerical computations

A Matlab implementation of the path-sum approach, presented in this work, has been made freely available [22]. Here we present an assessment of its computational performances compared with the standard piecewise-constant propagator approximation, based on the identity

U ptq " lim ∆tÑ0 N"t{∆t ź k"0 e ´iH pk∆tq∆t , (28)
and adaptive Runge-Kutta method, using ode45 in Matlab.

Exemplar waveform: chirped pulse

In this article we employ a general expression of βptq for a chirped pulse [START_REF] Foroozandeh | Spin dynamics during chirped pulses: applications to homonuclear decoupling and broadband excitation[END_REF], which can be expressed using 6 parameters: amplitude (ω 1 ), bandwidth (∆F), duration (τ p ), overall phase (φ 0 ), time offset (δ t ), and frequency offset (δ f ). The time envelope (amplitude profile) of a chirped pulse can be expressed using a super-Gaussian distribution

ω 1 ptq " ω 1,max exp " ´2n`2 ˆt ´δt τ p ˙n , ( 29
)
where n is a smoothing factor. The frequency sweep function can be written as

ωptq " 2π∆Fpt ´δt q τ p `2πδ f , (30)
which correspond to the phase

φptq " φ 0 `π∆Fpt ´δt q 2 τ p `2πδ f pt ´δt q. ( 31
)
The maximum amplitude of a chirped pulse (ω 1,max ) can be calculated using three parameters as

ω 1,max " d 2π∆FQ 0 τ p , (32)
where Q 0 is the adiabaticity factor at time τ p{2, where t P r0, τ p s and δ t " τ p{2. The relationship between Q 0 and pulse flip angle (α) can be written as [START_REF] Jeschke | Coherence transfer by passage pulses in electron paramagnetic resonance spectroscopy[END_REF] Q

0 " 2 π ln ˆ2 cos pαq `1 ˙. (33)
As the effective flip angle of a chirped pulse approaches 180 asymptotically as the ω 1 increases, for most practical purposes a value of Q 0 is chosen to satisfy adiabatic condition, while avoiding excessive pulse amplitudes.

Figure 4 shows comparisons between numerical and analytical results for single spin-1 {2 under a chirped pulse. Figure 5 shows the performance of several analytical approximations as in eq. ( 18), for a chirped pulse with realistic parameters, versus the true solution. Figure 6 shows the time-evolution of a spin-1 2 experiencing an adiabatic inversion on a Bloch sphere, obtained from the path-sum solution described in Section 3.1.

Numerical implementation of path-sums

Let I " r0, T s be the time interval over which the numerical solution is sought. Let tt i P Iu 0ďiďN´1 be the discrete times at which the solution is to be numerically evaluated. For simplicity, we take these time points to be equally spaced by ∆t. For a smooth function f pt 1 , tq over I 2 , we define the triangular matrix F with entries

F i, j :" # f pt i , t j q, for i ě j 0, otherwise.
With discretized time, the Volterra composition ‹ of eq. ( 13) turns into an ordinary matrix product and as predicted by Equations ( 14) and ( 16) (dashed red line), for a monopartite system with Ω " 2π 7kHz under a chirped pulse, described in Section 5.1, with ω 1 " 2π ˆ1545 rad/s, ∆F " 30kHz, τ p " 10ms, φ 0 " 0, δ t " 5ms, δ f " 0, and n " 30. Two curves show perfect overlap and hence are indistinguishable. Inset shows the normalised temporal amplitude profile 2|βptq|{ω 1 .

p f ‹ gqpt i , t j q " ż t i t j f pt i , τqgpτ, t j qdτ § đ ÿ t j ďt k ďt i f pt i , t k qgpt k , t j q∆t " pF .G q i, j ∆t, (34)
U 22 (1) (t)
U 22 (30) (t) where G is the triangular matrix built from g. Most importantly, this observation extends to ‹-resolvents

`1‹ ´f ˘‹´1 pt i , t j q " lim ∆tÑ0 1 ∆t `I ´∆t F ˘´1 i, j .
Therefore, in practical numerical computations with ∆t ! 1 the ordinary inverse of the triangular matrix I ´∆t F approximates the ‹-resolvent of f . Furthermore, matrix I ´∆t F is necessarily well conditioned since its diagonal entries 1 ∆t F i,i can be made arbitrarily close to 1 by choosing ∆t small enough. Consequently, numerical evaluation of any path-sum requires only multiplying and inverting well-conditioned triangular matrices [START_REF] Birkandan | Computations of general Heun functions from their integral series representations[END_REF][START_REF] Giscard | Lanczos-like method for the time-ordered exponential[END_REF] We improve upon this strategy by noting that eq. ( 34) corresponds to using the rectangular rule of integration. Using the trapezoidal or averaged Simpson rule [START_REF] Kalambet | Comparison of integration rules in the case of very narrow chromatographic peaks[END_REF] instead leads to much more accurate results. Standard numerical analysis indicates that a code using trapezoidal quadrature on N time points and time step ∆t " T {N should have an accuracy scaling as Op∆t ´2q and a computational cost of OpN 2 q. Similarly, a code relying on the average Simpson rule of integration should have an accuracy scaling as Op∆t 3 q for a computational cost of OpN 2 q. In contrast, the piecewise-constant propagator approximation of eq. ( 28) is expected to have accuracy Op∆tq and linear computational cost OpNq.

Seeking more flexibility in the calculations, the path-sum codes allow for subdivisions of the time interval I into N I smaller intervals. Over each of the subintervals, the evolution operator is evaluated from its path-sum formulation on N p discrete time points. The evolution operator at the end of the interval is passed as a seed to the next interval, over which it is once again evaluated from its path-sum. This offers an hybrid approach between piecewise-constant propagator approximation (PCPA) (N p " 1, N I " 1) and pure path-sum pN p " 1, N I " 1). Overall, this hybrid approach evaluates the solution at a total of N " N p N I time points. Tuning N p and N I independently allows the user to trade accuracy for speed and vice-versa. For example, keeping N p moderate while choosing N I " 1 allows for faster evaluations than N p " 1, N I " 1 and is thus better suited should the user require the solution at a very large number of time points. At the opposite, more accurate results will be obtained by making N p large while N I can be as low as 1. In practise, we found that in most-though not all-situations a pure path-sum approach (N p " 1, N I " 1) offers the best trade-off of speed versus accuracy.

The codes designed to produce the entire discretized timeevolution of the density matrix taking the spin system and the waveform parameters as inputs. Figure 7 shows the output density matrix for a bipartite (4 ˆ4) system undergoing an adiabatic inversion using a chirped pulse.

Accuracy evaluation

Let ρ M ptq be the density matrix as evaluated by a method M (e.g. path-sum or PCPA) and let ρ r be the reference density matrix as evaluated within machine precision (relative and absolute tolerances set to 10 ´13 for each entry) by the standard solver ode45. We evaluate the accuracy of ρ M ptq by evaluating the deviation from 1 of its normalised Frobenius scalar product with ρ r Example of an output produced by the path-sum code using Simpson quadrature in the case of a bipartite system with relative error E M " 3.28 10 ´9 using eq. ( 35). Each element of the figure represents the real part of the density matrix elements ρ i, j ptq as a function of time t (in µs). Spin system parameters: Ω 1 " 2π ˆ700 rad/s, Ω 2 " 2π ˆ600 rad/s, J " 150 Hz; waveform parameters: τ p " 1ms, ∆F " 50kHz, φ 0 " 0, ω 1 " 2π ˆ6.31 krad/s, δ t " 0.5ms, δ f " 0 and n " 20.

E M :" 1 T ż T 0 1 ´Tr `ρ: M pτqρ r pτq ȃ}ρ M pτq} F }ρ r pτq} F dτ, (35)
which is the relative error on ρ M with respect to the reference solution. Here }A} F :" TrpA : Aq designates the Frobenius norm of matrix A. As constructed above, the relative error E M evaluates to 0 if ρ M ptq " ρptq at all times.

Performances

Throughout this section we consider a βptq for a chirped pulse as described in Section 5.1. All computational tests take place on a Dell Laptop running Ubuntu 18.04 equipped with Intel Core i7-8665U CPU @ 1.90GHz 8 running Matlab R2019a. In Tables 1 to 3 we show the time and number of points N needed for each method to reach a target relative error E M of 10 ´3, 10 ´6 and 10 ´8. In the cases of Path-Sum (PS) codes, we have always set N p " N and N I " 1 to facilitate comparisons.

These tests show that for relative errors E M on the order of 10 ´3 to 10 ´7-largely sufficient to simulate actual experimentspath-sum codes outperform the piecewise-constant propagator approximation (see in particular Figure 8) and even ode45. The latter observation is surprising given that: i) hear we are dealing with small systems for which ode45 is highly optimized; ii) ode45 is allowed to choose the number and position of its computation points while path-sum is constrained by our implementation to work on equally spaced points-a known disadvantage; iii) ode45's computation time degrades rapidly when asked to produce the solution at a large number of evaluation points.

Further improvements to the current path-sum implementation can be made by using optimized, instead of linearly spaced, time points or working with an orthogonal polynomial basis. Additionally, C or Python implementations can offer significant speed-up compared to the present Matlab one, as evidenced by path-sum based codes for computing Heun functions [START_REF] Birkandan | Computations of general Heun functions from their integral series representations[END_REF]26]

Numerical outlook for large systems

As the size of the Hamiltonian increases, we expect the difference of computation times to grow in favor of codes evaluating the path-sum solutions, in particular, for codes exploiting path-sum's scale invariance property [START_REF] Giscard | Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR[END_REF] or using Lanczos path-sum [START_REF] Giscard | Lanczos-like method for the time-ordered exponential[END_REF][START_REF] Giscard | Tridiagonalization of systems of coupled linear differential equations with variable coefficients by a Lanczos-like method[END_REF][START_REF] Giscard | Lanczos-Like algorithm for the time-ordered exponential: The ˚-inverse problem[END_REF]. Lanczos path-sum is a type of preconditioning procedure for path-sum that finds a tridiagonal time-dependent matrix T ptq whose line 1 column 1 entry of the time-ordered exponential is the same as v T .U ptq.w, with any pair of vectors, v and w, i.e. T is determined such that T exp ˆż t 0 T pτqdτ ˙1,1 " v T .U ptq.w.

While the above equality is exactly satisfied for a tridiagonal matrix T of the same size as U , an excellent approximation of v T .U ptq.w may be achieved using a truncation of T that is much smaller than U . This procedure, standard in timeindependent Lanczos methods, partially alleviates the problem posed by the exponential expansion of the quantum space with the number of spins the Hamiltonian's sparsity and can effectively reduce the size of the matrices involved while maintaining numerical accuracy by truncating T ptq [START_REF] Giscard | Lanczos-like method for the time-ordered exponential[END_REF]. Once a suitable truncation of T is determined, path-sum is exploited to compute its time ordered exponential. Because T and its truncations are all tridiagonal the path-sums are finite scalar-valued ‹-continued fractions with a single branch similar in structure to those presented in this work.

Conclusion

In this work we used the path-sum formalism and showed the wide applicability of this approach to solve, both analytically and numerically, the spin dynamics of quantum spin systems driven by time dependent forces. Analytically speaking, the method offers closed form expressions in terms of ‹-resolvents, which may nonetheless themselves be transcendent mathemati-cal functions. However, these ‹-resolvents are always available from their analytical unconditionally converging Neumann series expansions or numerically, to any desired accuracy. Using the time-discretization of Volterra compositions, numerical implementations of path-sums require only multiplying and inverting well conditioned triangular matrices. Via a diverse set of examples, we demonstrated that the resulting numerical path-sum integrator consistently outperforms piecewiseconstant propagator approximation and ODE integration methods, which furthermore do not provide any analytical insights in contrast with the path-sum approach.

The proposed method has the potential to find applications in a variety of areas where having access to exact solutions of the time evolution of quantum spin systems is beneficial, including geometric [START_REF] Khaneja | Sub-riemannian geometry and time optimal control of three spin systems: Quantum gates and coherence transfer[END_REF][START_REF] Khaneja | Subriemannian geodesics and optimal control of spin systems[END_REF][START_REF] Jurdjevic | Optimal control on lie groups and integrable hamiltonian systems[END_REF][START_REF] Bonnard | A review of geometric optimal control for quantum systems in nuclear magnetic resonance[END_REF][START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF] and adiabatic optimal control [START_REF] Kato | On the adiabatic theorem of quantum mechanics[END_REF][START_REF] Brif | Exploring adiabatic quantum trajectories via optimal control[END_REF][START_REF] Meister | Optimal control theory with arbitrary superpositions of waveforms[END_REF][START_REF] Augier | Adiabatic ensemble control of a continuum of quantum systems[END_REF] methods, in addition to optimal control of quantum spin systems [START_REF] Meister | Optimal control theory with arbitrary superpositions of waveforms[END_REF].

p8q

  22 ptq. Nonetheless, the Neumann approximations U pmq 22 ptq are useful in that they provide analytical closed-form expressions for any desired accuracy.

  2πJ i j I piq ¨Ip jq , and I piq ¨Ip jq " I

Figure 4 :

 4 Figure 4: Time evolution of U 22 ptq as obtained numerically (solid black line) and as predicted by Equations (14) and (16) (dashed red line), for a monopartite system with Ω " 2π 7kHz under a chirped pulse, described in Section 5.1, with ω 1 " 2π ˆ1545 rad/s, ∆F " 30kHz, τ p " 10ms, φ 0 " 0, δ t " 5ms, δ f " 0, and n " 30. Two curves show perfect overlap and hence are indistinguishable. Inset shows the normalised temporal amplitude profile 2|βptq|{ω 1 .

Figure 5 :

 5 Figure 5: Exact time evolution of U 22 ptq (solid black line, corresponding to U p8q 22 ptq) and approximated evolutions, as in eq. (18), predicted by the analytical closed-forms Neumann approximations: U p1q 22 ptq (dotted red line), U p30q 22 ptq (dashed blue line). Inset shows the normalised temporal amplitude profile 2|βptq|{ω 1 . Spin and pulse parameters are identical to those given in the caption of Figure 4.

Figure 6 :

 6 Figure 6: Time evolution |ψptqy " U : .U ptq.U.ψ 0 on the Bloch sphere as dictated by the path-sum solution described Section 3.1. Spin and pulse parameters are identical to those given in the caption of Figure 4.

Figure 7 :

 7 Figure7: Example of an output produced by the path-sum code using Simpson quadrature in the case of a bipartite system with relative error E M " 3.28 10 ´9 using eq. (35). Each element of the figure represents the real part of the density matrix elements ρ i, j ptq as a function of time t (in µs). Spin system parameters: Ω 1 " 2π ˆ700 rad/s, Ω 2 " 2π ˆ600 rad/s, J " 150 Hz; waveform parameters: τ p " 1ms, ∆F " 50kHz, φ 0 " 0, ω 1 " 2π ˆ6.31 krad/s, δ t " 0.5ms, δ f " 0 and n " 20.

Figure 8 :

 8 Figure 8: Relative error E M versus computation time for PCPA (red disks), path-sum trapezoidal ('PS-T', blue squares) and path-sum Simpson ('PS-S', black diamonds) in the case of the bipartite system with the parameters of Figure 7.

  [START_REF] Hioe | Solution of bloch equations involving amplitude and frequency modulations[END_REF] ptq "

			which leads to
							U 12 ptq "	´i ?	2e iΩt	ż t	e ´iΩτ βpτqU 22 pτqdτ,
							0
			and G 32 ptq " Ğ G 12 ptq. Here U 22 and G 22 are given by Equa-
			tions (14) and (16) respectively. Similarly we get
			G 21 ptq "
							Simple path 2Ð1
			hkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkikkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkj p1 ‹ ´i ? 2β ‹ p1 ‹ `iΩq ‹ i ? 2 βq ‹´1 loooooooooooooooooooooomoooooooooooooooooooooon ‹ ´i ? 2 β loomoon ‹ lo omo on G 11 ,
						All walks 2Ð2 avoiding 1	Edge 2Ð1	All walks 1Ð1
			"	´i ?	2	ż t ´1‹ `2e ´iΩt 1 βpt 1 q Ě B t 1 ,t pΩ{2πq ¯‹´1 pt, τq βpτqU 11 pτqdτ,
							0
			and U 21 "	ş t 0 G 21 pτqdτ. By symmetry, U 23 ptq " Ğ U 21 ptq. We
			conclude with
							Simple path 3Ð2Ð1
				G 31 ptq "	hkkkkkkkkkkkkkkkkkkkkkikkkkkkkkkkkkkkkkkkkkkj p1 ‹ `iΩq ‹´1 loooooomoooooon ‹ i ? 2 β lo omo on ‹ lo omo on G 21 ,
							Loop 3Ð3	Edge 3Ð2	All walks 2Ð1
			which yields
							U 31 ptq " i	?	2e ´iΩt	ż t	e iΩτ	s βpτq U 21 pτqdτ.
							0
			and, by symmetry, U 13 ptq " Ğ U 31 ptq.
			Just as for the diagonal terms, all of the path-sum results for
			the off-diagonal terms yield unconditionally convergent Neu-
			mann expansions, e.g.
			G 21 ptq "	´i ?	2 βptqU 11 ptq
							`2i	?	2 e ´iΩt 1 βpt 1 q	ż t	Ě B t,τ pΩ{2πq βpτqU 11 pτqdτ `¨¨4
							0
			. Solutions in SUp2 N q
			4.1. Liouvillevon Neumann equation
			The solution of the Liouvillevon Neumann equation in the
			Hilbert space can be written using a time-dependent unitary
			propagator as
							ρptq " U ptqρp0qU ptq : ,
		Simple path 1Ð2
	hkkkkkkkkkkkkkkkkkkkkkkkkkkkkikkkkkkkkkkkkkkkkkkkkkkkkkkkkj p1 ‹ ´p´iqp´Ωq loooomoooon q ‹´1 ‹ p´i ? 2βq looomooon lo omo on ‹ G 22 ,
	Loop 1Ð1	Edge 1Ð2	All walks 2Ð2

  1 {2 system with offsets Ω 1 and Ω 2 and a coupling constant J can be expressed as

	H ptq "	¨h11 βptq βptq βptq h 22 πJ βptq 0 ‹ ‹ βptq πJ h 33 βptq '
		0	βptq βptq h 44

  ż t 1 t e ´iH II pt 1 ´τq .

				ˆβpτq βpτq βpτq βpτq	˙dτ	¸‹´1	,
	and
	ˆG2,1 , G 3,1	˙"
	G 2,4 , G 3,4
	ż t 0	e ´iH II pt´τq .	βpτq βpτq ˆβpτq βpτq ˙.ˆU U 4,1 pτq U 4 pτq 1 pτq U 1,4 pτq ˙dτ.

  H II " ¨h22 πJ 23 πJ 13 πJ 23 h 33 πJ 12 πJ 13 πJ 12 h 55 ¨h44 πJ 12 πJ 13 πJ 12 h 66 πJ 23 πJ 13 πJ 23 h 77

	',
	H III "

Table 1 :

 1 . Comparisons of computational performances for a monopartite system with Ω " 2π ˆ1 krad/s in SU(2) representation driven by a chirped pulse with parameters: τ p " 1ms, ∆F " 100 kHz, φ 0 " 0, ω 1 " 2π ˆ8.92 krad/s, δ t " 0.5ms, δ f " 0, n " 30 over the time interval up to T " 1ms.

			Method M	E M	N	Time (s)
			ode45	10 ´3	´0.092
			PCPA	10 ´3	500	0.087
		PS Trapezoidal 10 ´3	140	0.01
			PS Simpson	10 ´3	121	0.007
			ode45	10 ´6	´0.113
			PCPA	10 ´6	17000	1.596
		PS Trapezoidal 10 ´6	2000	0.115
			PS Simpson	10 ´6	300	0.012
			ode45	10 ´8	´0.121
			PCPA	10 ´8 450000	42.86
		PS Trapezoidal 10 ´8	15000	2.69
			PS Simpson	10 ´8	7000	2.49
		0.001
				PCPA
	Relative error	10 -5	PS-T
		10 -7	PS-S
		10 -9
			0.05 0.10		0.50	1	5	10
				Computation time (s)

  . The ‹-Lanczos approach heavily exploits

	Method M	E M	N	Time (s)
	ode45	10 ´3	´0.092
	PCPA	10 ´3	510	0.011
	PS Trapezoidal 10 ´3	127	0.022
	PS Simpson	10 ´3	120	0.017
	ode45	10 ´6	´0.095
	PCPA	10 ´6	16000	0.040
	PS Trapezoidal 10 ´6	700	0.251
	PS Simpson	10 ´6	350	0.035
	ode45	10 ´8	´0.102
	PCPA	10 ´8 160000	0.255
	PS Trapezoidal 10 ´8	2000	4.96
	PS Simpson	10 ´8	3000	6.18

Table 2 :

 2 Comparisons of computational performances for a monopartite system in SO(3) representation driven by a chirped pulse with parameters given in the caption of Table1. Here path-sum codes solve the whole 3 ˆ3 system of Bloch equations irrespective of the initial state, ψp0q. In contrast, the PCPA-based algorithm evolves ψptq from one given ψp0q, an easier task.

	Method M	E M	N	Time (s)
	ode45	10 ´3	´0.172
	PCPA	10 ´3	330	0.075
	PS Trapezoidal 10 ´3	85	0.030
	PS Simpson	10 ´3	77	0.026
	ode45	10 ´6	´0.221
	PCPA	10 ´6	10500	1.19
	PS Trapezoidal 10 ´6	500	0.206
	PS Simpson	10 ´6	200	0.046
	ode45	10 ´8	´0.269
	PCPA	10 ´8 110000	12.47
	PS Trapezoidal 10 ´8	2000	8.89
	PS Simpson	10 ´8	1100	1.69

Table 3 :

 3 Comparisons of computational performances in the bipartite cases with parameters given in the caption Figure7.

Acknowledgements

PLG acknowledges funding from the Agence Nationale de la Recherche, grants ANR-19-CE40-0006 and ANR-20-CE29-0007. MF is grateful to the Royal Society for a University Research Fellowship and a University Research Fellow Enhancement Award (grant numbers URF\R1\180233 and RGF\EA\181018). We thank E. Baligacs for an improvement to the Matlab codes and for suggesting the Frobenius scalar product to evaluate accuracy.


A S Vasenko

S Kawabata

A A Golubov

M Yu Kupriyanov

Claudine Lacroix

F S Bergeret

Frank W J Hekking

Current-voltage characteristics of tunnel Josephson junctions with a ferromagnetic interlayer

Keywords: numbers: 74.45.+c, 74.50.+r, 74.78.Fk, 75.30.Et

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The possibility to switch the ground-state of a Josephson junction from a 0 to a π phase state and the possible application of such junctions in quantum information led to a renewal interest in the study of the so called π Josephson junctions. The existence of such a transition was predicted more than thirty years ago, 1 however due to technological requirements only recently it was observed. The realization of π Josephson junctions was achieved in superconductor/ ferromagnet/ superconductor (SFS) junctions. [2][3][4][5][6][7][8][9][10][11][12][13][14][START_REF] Weides | Josephson Tunnel Junctions with Ferromagnetic Interlayer[END_REF][START_REF] Pfeiffer | [END_REF][17][18][19][20][21] Microscopically, S/F hybrid structures are characterized by an unusual proximity effect, with a damped oscillatory behavior of the superconducting correlations in the F layer (for a review see Refs. 22-24 and references therein). This unusual proximity effect in S/F layered structures leads to a number of striking phenomena like the nonmonotonic dependence of their critical temperature and the appearance of oscillations of critical current in SFS Josephson junctions as a function of the F layer thickness. 2,9 In particular the change of sign of the critical current corresponds to the so-called 0 -π transition.

On the other hand, SFS junctions, as any metallic junction exhibit very small resistances and therefore are not quite suitable for those applications, for which active Josephson junctions are required. This problem can be solved by adding an additional insulating (I) layer to increase the resistance. SIFS junctions represent an interesting case for practical use of π Josephson junctions. For instance, a SIFS structure offers the freedom to tune the critical current density over a wide range and at the same time realize high values of the product of the junction critical current I c and its normal state resistance R n . [13][14][START_REF] Weides | Josephson Tunnel Junctions with Ferromagnetic Interlayer[END_REF] In addition, Nb based tunnel junctions are usually underdamped, which is desired for many applications. Due to these advantages, SIFS π junctions have been proposed as potential elements in superconducting classical and quantum logic circuits. 25,26 For instance, SIFS junctions can be used as complementary elements (π -shifters) in RSFQ circuits (see Ref. 27 and references therein). Finally, SIFS structures have been proposed for the realization of so called ϕ-junctions with a ϕ drop in the ground state, where 0 < ϕ < π. 28 The properties of SIFS junctions have been intensively studied both experimentally [9][10][11][12][13][14][START_REF] Weides | Josephson Tunnel Junctions with Ferromagnetic Interlayer[END_REF][START_REF] Pfeiffer | [END_REF][17][18] and theoretically. [28][29][30][31] However, properties of the quasiparticle current have received relatively little attention so far, although they can be very important for the description of SIFS junctions as possible elements of superconducting logic circuits.

The purpose of this work is to provide a quantitative model describing the behavior of quasiparticle current in SIFS junctions as a function of parameters characterizing material properties of the ferromagnetic interlayer. We also focus our study on the properties of the density of states (DOS) in S/F bilayers and discuss the oscillations of the DOS at Fermi energy. Finally, we calculate the macroscopic quantum tunneling (MQT) escape rate for current-biased SIFS junctions by taking into account the dissipative correction due to the quasiparticle tunneling. Based on this we conclude that the influence of the quasiparticle dissipation on the macroscopic quantum dynamics of SIFS junctions is small, which is an advan-FIG. 1: (Color online) Geometry of the considered system. The thickness of the ferromagnetic interlayer is d f . The transparency of the left (right) S/F interface is characterized by the coefficient γ B1(B2) . The left interface is an insulating barrier, γ B1 ≫ 1 (shown by a black line), while the right interface is transparent, γ B2 ≪ 1 (shown by a grey line).

tage of SIFS junctions for quantum logic (qubit) applications.

The paper is organized as follows. In the next section we formulate the theoretical model and the basic equations. In Sec. III we solve the nonlinear Usadel equations numerically for arbitrary length of the ferromagnetic layer, and calculate the DOS in the F layer. We compare these results with an analytical expression for DOS in case of a long SIFS junction, i.e. when the thickness d f of the ferromagnetic layer is much larger than the decay length of the characteristic superconducting correlations in the ferromagnet ξ f 1 . We also discuss the oscillations of the DOS at the Fermi energy. In Sec. IV we present the current-voltage characteristics of SIFS junctions for different parameters of the ferromagnetic interlayer. In Sec. V, we use these data to calculate the MQT escape rate for current-biased SIFS junctions by taking into account the dissipative effect of the quasiparticle tunneling. Finally we summarize the results in Sec. VI.

II. MODEL AND BASIC EQUATIONS

We consider a SIFS junction such as the one depicted in Fig. 1. It consists of a ferromagnetic layer of thickness d f and two thick superconducting electrodes along the x direction. The left and right superconductor/ ferromagnet interfaces are characterized by the dimensionless parameters γ B1 and γ B2 , 32,33 respectively, where γ B1,B2 = R B1,B2 σ n /ξ n , R B1,B2 are the resistances of the left and right S/F interfaces, respectively, σ n is the conductivity of the F layer, ξ n = D f /2πT c , D f is the diffusion coefficient in the ferromagnetic metal and T c is the critical temperature of the superconductor (we assume h = k B = 1, except for Sec. V). We also assume that the S/F interfaces are not magnetically active. We will consider the diffusive limit, in which the elastic scattering length ℓ is much smaller than the decay characteristic length ξ f 1 = min ξ f 1↑ , ξ f 1↓ (for the definitions of ξ f 1↑(↓) see Eqs. (16) be- low).

We assume that the tunneling barrier is located at the left S/F interface, while the right interface is perfectly transparent, this means that γ B1 ≫ 1 while γ B2 ≪ 1. In this case the left S layer and the right F/S bilayer in Fig. 1 are decoupled and we can calculate the quasiparticle current through a SIFS junction using the standard tunneling formula 34

I = 1 eR ∞ -∞ dEN s (E -eV )N f (E) [ f (E -eV ) -f (E)], (1)
where

N s (E) = |E|Θ(|E| -∆)/ √ E 2 -∆ 2 is the BCS DOS, Θ(x) is the Heaviside step function, N f (E) is the DOS in the ferromagnetic interlayer at x = -d f /2, f (E) = [1 + exp(E/T )] -1 is the Fermi function, and R ≡ R B1 . Both N s (E)
and N f (E) are normalized to their values in the normal state. In particular at zero temperature, T = 0, the current acquires the form,

I = Θ(eV -∆) 1 eR eV -∆ 0 dEN s (E -eV )N f (E) . (2)
To obtain N f (E) we notice that since γ B1 ≫ 1, the left su- perconducting lead does not influence the DOS in the ferromagnetic interlayer (to zero order in the barrier transparency). This reduces the problem to the following: we need to find the DOS of a single F/S bilayer, which can be done by solving the Usadel equations in the ferromagnetic layer.

Using the θ -parameterizations of the normal and anomalous Green functions, G = cos θ , F = sin θ , we can write the Usadel equations in the F layer as [35][36][37]

D f 2 ∂ 2 θ f ↑(↓) ∂ x 2 = ω ± ih + 1 τ z cos θ f ↑(↓) sin θ f ↑(↓) + 1 τ x sin θ f ↑ + θ f ↓ ± 1 τ so sin θ f ↑ -θ f ↓ , (3)
where the positive and negative signs correspond to the spin up ↑ and spin down ↓ states respectively. In this notation the spin up state corresponds to the anomalous Green function F ↑ ∼ ψ ↑ ψ ↓ while the spin down state corresponds to F ↓ ∼ ψ ↓ ψ ↑ , where ψ ↑(↓) are the electron fermionic opera- tors. The ω = 2πT (n + 1 2 ) are the Matsubara frequencies, and h is the exchange field in the ferromagnet. The scattering times are labelled here as τ z , τ x and τ so , where τ z(x) cor- responds to the magnetic scattering parallel (perpendicular) to the quantization axis, and τ so is the spin-orbit scattering time. [38][39][40][41] We consider here ferromagnets with a strong uniaxial anisotropy, in which case the magnetic scattering does not couple the spin up and spin down electron populations, i.e. the perpendicular fluctuations of the exchange field are suppressed (τ -1

x ∼ 0). Therefore, we will neglect τ x in our con- sideration and denote τ z as a magnetic scattering time τ m . We will also consider ferromagnets with weak spin-orbit interactions and henceforth also neglect the spin-orbit scattering time τ so . In this case the Usadel equations in the ferromagnetic layer for different spin projections are not coupled any more and can be written as, 4) while in the S layer the Usadel equations take the form

D f 2 ∂ 2 θ f ↑(↓) ∂ x 2 = ω ± ih + cos θ f ↑(↓) τ m sin θ f ↑(↓) , (
D s 2 ∂ 2 θ s ∂ x 2 = ω sin θ s -∆(x) cos θ s . (5)
Here D s is the diffusion coefficient in the superconductor and ∆(x) is the superconducting pair potential. Notice that in the latter equation we have omitted the subscripts '↑ (↓)' because both equations are identical in the superconductor. Eqs. ( 4), ( 5) should be complemented by the selfconsistency equation for the superconducting order parameter ∆, 6) and by the boundary conditions at the outer boundary of the ferromagnet

∆(x) ln T c T = πT ∑ ω>0 2∆(x) ω -sin θ s↑ -sin θ s↓ , (
∂ θ f ∂ x -d f /2 = 0, (7)
and at the F/S interface, 32

ξ n γ ∂ θ f ∂ x d f /2 = ξ s ∂ θ s ∂ x d f /2 , ( 8a
)
ξ n γ B2 ∂ θ f ∂ x d f /2 = sin θ s -θ f d f /2 , (8b)
where γ = ξ s σ n /ξ n σ s , σ s is the conductivity of the S layer and ξ s = D s /2πT c . The parameter γ determines the strength of suppression of superconductivity in the right S lead near the interface compared to the bulk: no suppression occurs for γ = 0, while strong suppression takes place for γ ≫ 1. In our numerical calculations we will assume small γ ≪ 1. Notice that the interface parameters do not depend on the spin direction. In other words we are not considering spin-active interfaces. In the case of spin-active barriers, one should use the boundary conditions introduced in Refs. 42-44, rather than the standard Kupriyanov-Lukichev boundary conditions [Eqs. (8)].

To complete the boundary problem we also set a boundary condition at

x = ∞, θ s (∞) = arctan ∆ ω , (9)
where the Green functions acquire the well known bulk BCS form. Eqs. ( 4)-( 9) represent a closed set of equations that should be solved self-consistently. As it will be discussed in the next section, the knowdlege of the Green function will allow us to compute the DOS at the outer F boundary.

III. DENSITY OF STATES IN THE F/S BILAYER

The DOS N f (E) normalized to the DOS in the normal state, can be written as

N f (E) = N f ↑ (E) + N f ↓ (E) /2, (10)
where N f ↑(↓) (E) are the spin resolved DOS written in terms of spectral angle θ ,

N f ↑(↓) (E) = Re cos θ f ↑(↓) (iω → E + i0) . (11)
To obtain N f , we use a self-consistent two-step iterative procedure 37,[45][46][47] . In the first step we calculate the pair potential coordinate dependence ∆(x) using the self-consistency equation in the S layer, Eq. ( 6). Then, by proceeding to the analytical continuation in Eqs. ( 4), ( 5) over the quasiparticle energy iω → E + i0 and using the ∆(x) dependence obtained in the previous step, we find the Green functions by repeating the iterations until convergency is reached.

Before showing the numerical results we consider an analytic limiting case. If the F layer is thick enough (d f ≫ ξ f 1 ) and γ = 0 in Eq. (8a), the DOS at the free boundary of the ferromagnet can be written as 29,48

N f ↑(↓) (E) = Re[cos θ b↑(↓) ] ≈ 1 - 1 2 Re θ 2 b↑(↓) . (12)
Here θ b↑(↓) is the value of θ f at x = -d f /2, given by

θ b↑(↓) = 8F(E) (1 -η 2 )F 2 (E) + 1 + 1 exp -p d f ξ f , (13)
where ξ f = D f /h. In Eq. ( 13) we use the following notations,

p ↑(↓) = 2/h -iE R ± ih + 1/τ m , (14a) η 2 ↑(↓) = (1/τ m )(-iE R ± ih + 1/τ m ) -1 , (14b)
F(E) = ∆ -iE R + ∆ 2 -E 2 R , E R = E + i0. (14c)
Here, we again adopt the convention that a positive (negative) sign in front of h corresponds to the spin up state ↑ (spin down state ↓). Hereafter we will write spin labels ↑ (↓) explicitly only when needed. From Eqs. ( 12)-( 13) we obtain for the full DOS the following expression in the limit

d f ≫ ξ f 1 , N f ≈ 1 -Re ∑ ↑, ↓ 16F 2 (E) exp -p 2d f ξ f ( (1 -η 2 )F 2 (E) + 1 + 1) 2 . ( 15
)
At this point, we define the characteristic decay and oscillation lengths ξ f 1,2↑(↓) as

p ↑(↓) /ξ f = 1/ξ f 1↑(↓) + i sgn(h ∓ E)/ξ f 2↑(↓) , (16a)
1 ξ f 1↑(↓) = 1 ξ f E ∓ h h 2 + 1 h 2 τ 2 m + 1 hτ m , ( 16b
) 1 ξ f 2↑(↓) = 1 ξ f E ∓ h h 2 + 1 h 2 τ 2 m - 1 hτ m . (16c)
In the absence of magnetic scattering both lengths coincide and are equal to ξ f h/|E ∓ h| for different spin orientations. One can rewrite Eq. ( 15) in the following form,

N f ≈ 1 -∑ ↑, ↓ exp - 2d f ξ f 1 A sin χ + 2d f ξ f 2 + B cos χ + 2d f ξ f 2 , (17) 0 4
1.00

1.05 N f (E) E/ 0 4 0.8 1.0 N f (E) E/ 0 1 2 0.8 1.6 N f (E) E/ 0.5 1.0 1.5 0 2 N f (E) E/ c) d) b) a)
d f /ξ n = 0.5, (b): d f /ξ n = 1, (c): d f /ξ n = 2, (d): d f /ξ n = 3.
The approximate analytical solution, Eq. ( 15), is shown by dashed red lines.

where the coefficients A, B and χ can be obtained by expansion of the real part in Eq. ( 15); only two of them are independent. This form explicitly shows the damped oscillatory behavior of superconducting correlations in the F layer. The lengths ξ f 1,2 are also the lengths of decay and oscilla- tions of the critical current in SIFS junctions (see Eqs. (26) in Ref. 29). The period of the DOS oscillations is approximately twice smaller than the period of the critical current oscillations and the exponential decay is approximately twice faster than the decay of the critical current. 29 Now we turn to the exact numerical solution. The obtained energy dependencies of the DOS at the free F boundary of the F/S bilayer are presented in Figs. 2,3,4. The exchange field is chosen such that h > ∆, which corresponds to the experimental situation.

Figure 2 shows the DOS energy dependence for different d f in the absence of magnetic scattering. At small d f we observe the DOS double peak due to the Zeeman splitting of the BCS peak at E = ∆. Most probably in the experiments, the BCS Zeeman splitted peak as presented in Fig. 2 (a) will be seen as a single peak due to many-body interaction effects, which introduce a finite lifetime (damping) of the quasiparticles. We also observe that at small d f and relatively small exchange field h, full DOS turns to zero inside a minigap, which vanishes with the increase of d f .

The minigap also exists in the normal metal (N) DOS in the S/N bilayers. If the thickness d n of the normal metal is larger than the coherence length, the characteristic scale of the minigap is set by the Thouless energy, E T h = D n /d 2 n , where D n and is the diffusion coefficient of the normal metal. 46 In the F layer of the S/F bilayer, the exchange field h shifts the DOS for the two spin subbands in opposite directions, therefore the critical value h c of the exchange field at which the minigap in the spectrum closes can be roughly estimated as 49

h c ∼ E T h , E T h = D f /d 2 f . ( 18
)
This equation shows the qualitative tendency: for smaller d f a higher h is needed to close the minigap [see also Fig. 3 (a)].

The estimation, Eq. ( 18), is only valid in the absence of magnetic scattering, since τ m also influences the minigap 50 [see also Fig. 4 (a)].

In Fig. 2 we also observe that after the minigap closes, the DOS at the Fermi energy N f (0) rapidly increases to values larger than unity with further increase of d f ; then it oscillates around unity while its absolute value exponentially approaches unity [see also Fig. 5]. This is the well-known damped oscillatory behavior of the DOS in F/S bilayers. Experimental evidence for such behavior was provided by Kontos et al. 51 In the case of long enough ferromagnetic layer we also observe the DOS peak at E = h, which was previously discussed in Ref. 52. A similar effect was also discussed in N/F/S structures, where it was shown that a zero energy peak appears in DOS if E T h = h. 53 We also show in Fig. 2 the analytical approximation, Eq. ( 15), which is in good agreement with the numerical result for thick enough ferromagnetic layers. In the numerically obtained curves the peak at E = h is smeared because of finite γ = 0.01 for the transparent F/S interface at x = d f /2.

In the absence of magnetic scattering we rewrite the analytical DOS expression, Eq. ( 15), for E ≥ ∆ in the following way,

N f (E) = 1 + ∑ ± 16∆ 2 cos 2d f ξ f |E±h| h (E + ε)( √ E + ε + √ 2ε) 2 e - 2d f ξ f |E±h| h , (19)
where

ε = √ E 2 -∆ 2 .
We can clearly see the exponential asymptotic of the peak at E = h from the Eq. (19). We should keep in mind that Eq. ( 19) is valid for large d f /ξ f , but nevertheless we may qualitatively understand why we do not see the peak at E = h for small ratio of d f /ξ f : if this factor is small the variation of the exponent exp{-2(d f /ξ f ) |E -h|/h} near the point E = h is also small. The peak is observable only for h of the order of a few ∆. For larger exchange fields the peak is very difficult to observe, since the energy dependent prefactor of the exponent in Eq. ( 19) decays as E -2 for E ≫ ∆.

Figure 3 shows the DOS energy dependence for different values of the exchange field h in the absence of magnetic scattering. For stronger exchange field the minigap closes at smaller d f , in qualitative correspondence with Eq. (18). From numerical calculations we obtain the following condition [see also Fig. 5 (c)], valid for τ -1 m ∼ 0, 54 In Fig. 3 we also observe the peak at E = h; at large enough exchange fields its amplitude can be much larger than the amplitude of the peak at E = ∆ (see Fig. 3 (d), blue dash-dotted curve). The existence of the DOS peak at E = h gives a possibility to measure the exchange field directly in experiment by measuring the F/S bilayer DOS in compounds with small magnetic scattering (since magnetic scattering is smearing the peak, see below). For example, in Ref. 55 were reported exchange fields for Pd 1-x Ni x with different Ni concentration, obtained by a fitting procedure. Considering Nb as a superconductor with ∆ =1.3 meV, we can estimate the exchange field in Pd 1-x Ni x : for 7% of Ni fitting gives h =2.8 meV, which is 2.2∆, and for 11.5% of Ni h =3.9 meV, which is 3∆. 55 It is interesting to use direct measurements of the DOS peak at E = h to check these fitting predictions of Ref. 55.

h c ≈ 0.77E Th ≈ 2.71∆ ξ n d f 2 . ( 20
)
Ferromagnetic metals with exchange fields of the order of few ∆ are crucially important for the fabrication of SIFS junctions, valid for superconducting logic applications. Presently used ferromagnets have h ≫ ∆, and therefore short oscillation length [see Eq. (16c)], which makes it difficult to control the F layer thickness. In already existing SIFS structures the roughness is often larger than desired precision of d f . 56 We hope that our results will trigger the experimental activity in finding ferromagnetic alloys with h of the order of few ∆.

Figure 4 shows the DOS energy dependence for different values of magnetic scattering time. Similarly to Fig. 3, for stronger magnetic scattering the minigap closes at smaller d f . Also the DOS peak at E = h, visible for long enough ferromagnetic layer, is smeared. The analytical solution (not shown), Eq. ( 15) also agrees quite well with the numerical results for d f ≫ ξ f 1 .

Although our results are obtained for weak ferromagnets, they can in certain cases be extended for ferromagnets with strong exchange fields, h ≫ ∆. In the absence of magnetic scattering the Usadel equation in energy representation, Eq. ( 4), can be rewritten as,

i 2 ∂ 2 θ f ↑(↓) ∂ y 2 = E h ∓ 1 sin θ f ↑(↓) , (21)
where y = x/ξ f is the dimensionless coordinate. In the case of h ≫ ∆, we can neglect the first term on the right hand side of Eq. ( 21) to obtain the subgap DOS. Thus, in that case the subgap structure scales with the length ξ f and for example the results presented in Fig. 2 for h = 4∆ also describe the DOS in the case of a high exchange field if one scales d f correspondingly. This procedure does not apply however for h = 2∆ (see Fig. 3), since in that case one cannot simply neglect the term E/h in Eq. ( 21).

To show explicitly the aforementioned DOS oscillations at the Fermi energy, 51 in Fig. 5 we plot the numerically calculated function

δ N f (d f , h, τ m ) = |1 -N f 0 |, N f 0 = N f (E = 0). (22)
Using Eqs. ( 12)-( 17) and ( 22) we get the analytical expression for the function δ N, valid for

d f ≫ ξ f 1 , δ N = 32 Re 1 ( 2 -η 2 0 + 1) 2 exp -p 0 2d f ξ f , (23)
where

p 0 = 2/h ih + 1/τ m , (24a) η 0 = (1/τ m )(ih + 1/τ m ) -1 . (24b)
At vanishing magnetic scattering we obtain,

δ N = 32 3 + 2 √ 2 cos 2d f ξ f exp - 2d f ξ f , (25)
in which case the characteristic lengths of decay and oscillations are equal to ξ f . The dependence of δ N f (d f ) on the ferromagnetic layer thickness d f at different values of exchange field and magnetic scattering time is presented in Fig. 5 (a) and (b). From Fig. 5 (a) we can see that with increasing exchange field h the minigap closes at smaller d f in agreement with Eq. ( 20), the period of the DOS oscillations at the Fermi energy decreases, and the damped exponential decay occurs faster. This is easy to see from Eq. ( 25), since in the absence of magnetic scattering δ N depends on h only as a function of ξ f . From Fig. 5 (b) we can see that with increasing α m = 1/τ m ∆ the period of the DOS oscillations on the contrary increases, although the minigap also closes at smaller d f and the damped exponential decay occurs faster. To understand this behavior we rewrite here the decay and oscillation lengths, Eq. ( 16), at the Fermi energy,

1 ξ f 1 = 1 D f h 2 + 1 τ 2 m + 1 τ m , (26a) 1 ξ f 2 = 1 D f h 2 + 1 τ 2 m - 1 τ m . ( 26b
)
We see from these equations that with increasing α m the length of decay ξ f 1 decreases, while the length of oscillations ξ f 2 increases. The dependencies of δ N f on exchange field and magnetic scattering time are presented in Fig. 5 (c) and (d), correspondingly. In Fig. 5 (c) we see oscillations of the DOS at the Fermi energy δ N f (h) around unity with increasing exchange field in the absence of magnetic scattering. In Fig. 5 (d) we show the function δ N(τ m ). It is interesting to note that its behavior can be both oscillatory and also monotonous. When the parameter α m increases starting from the minigap state (black solid curve) it is totally monotonous: increasing α m the minigap closes and the DOS starts to increase to unity, but never overshoots unity [we checked this up to α m = 80, see also Fig. 4 (a)]. If we start from the state where the minigap is already closed, we first observe oscillations, but then again a switch to monotonous behavior. For intermediate F layer thicknesses we see just one oscillation and then DOS monotonously approaches unity (we checked this up to α m = 80), while for thicker ferromagnets (d f /ξ n = 3) we observe two oscillations and then a monotonic behavior.

The dependencies δ N(h) and δ N(τ m ) can be important if in the experiment the material properties of the ferromagnetic interlayer, i.e. exchange field h and magnetic scattering time τ m can vary with some external parameter, for example tem- perature, magnetic field, etc.

Before turning to the calculation of the CVC we discuss briefly a recent experiment 57 in which a pronounced double peak in the DOS of Ni/Nb bilayers were reported. This double peak cannot be explained within our model based on the Zeeman splitting. The reason for the double peak in Ref. 57 remains controversial. In Ref. 44 it was numerically fitted by adding an extra parameter to the model, characterizing spinactive interfaces. However, this fit is far from being satisfactory. Nevertheless there is another feature of the DOS observed in Ref. 57 which can be explained within our model: by increasing d f the "normal" peak at E = ∆ [which is the BCS Zeeman split peak in Fig. 2(a)] is "inverted" [Fig. 2(b) and (c)] and becomes "normal" again [Fig. 2(d)] as d f is further increased. According to our model, at E = ∆ Eq. ( 19) reduces to the following expression,

N f (∆) = 1 + 16 ∑ ± cos 2d f ξ f h ± ∆ h e - 2d f ξ f h±∆ h . ( 27
)
This expression explains the inversion of the peak at E = ∆ as a function of d f . The peak is "normal" ("inverted") if the DOS at ∆ is larger (smaller) than unity. This variation is due to the sign of the cosine function in Eq. ( 27), which depends on the d f /ξ f ratio.

IV. CURRENT-VOLTAGE CHARACTERISTICS OF A SIFS JUNCTION

In this section we calculate the current-voltage characteristics (CVC) of a SIFS junction at low temperature, T = 0.1T c , using Eq. ( 1) and DOS N f (E, d f , h, τ m ) numerically obtained in the previous section.

Figure 6 shows the CVC of a SIFS junction in the absence of magnetic scattering. For comparison we also present the CVC of a SINS tunnel junction, i.e. a junction with a normal metal interlayer instead of a ferromagnet (h = 0). SINS structures were studied previously in Ref. 46. We observe several features of SIFS CVC which are the signatures of the proximity effect in the S/F bilayer.

For thin enough F layer we observe the "kink" on the CVC at eV ≈ 2∆ [Fig. 6 (a)], which corresponds to the case when the DOS N f exhibits a pronounced minigap. The corresponding DOS energy dependence (h/∆ = 4, d f /ξ n = 0.5) is shown in the inset. We can also see that for a certain range of parameters the CVC of a SIFS junction exhibit a nonmonotonic "wave" behavior. We can observe it for h/∆ = 4 (red short-dashed line) in Fig. 6 (b) and for h/∆ = 2 (blue dash-dotted line) in Fig. 6 (c). This behavior corresponds to the case when the DOS N f (E) minigap is already closed and the N f (0) at the Fermi energy is larger than unity. The corresponding DOS energy dependencies are presented in the insets of the plots 6 (b) and (c). At large enough d f and exchange fields the DOS N f (E) ≈ 1 and the current, Eq. ( 1), is given by the same equation as the current in the NIS tunnel junction,

0 2 0 2 IRe/ eV/ 1 2 0 2 IRe/ eV/ 0 2 0 2 IRe/ eV/ 1 2 0 2 IRe/ eV/ 0 1 0 1 N f (E) E/ 0 1 0 1 N f (E) E/ d) c) b) a)
I = 1 eR ∞ -∞ dEN s (E) [ f (E -eV ) -f (E)] . (28)
At small temperature T ≪ T c , this equation is well approximated by taking T = 0,

I = Θ(eV -∆) 1 eR eV ∆ E dE √ E 2 -∆ 2 = Θ(eV -∆) (eV ) 2 -∆ 2 eR , (29)
The red short-dashed line in Fig. 6 (d) (d f /ξ n = 3, h/∆ = 4, no magnetic scattering) almost coincides with this result except for the small region eV ≈ ∆, since for our numerically calculated curves we fix temperature T = 0.1T c ≪ T c . Figure 7 shows current-voltage characteristics of a SIFS junction in case of finite magnetic scattering. Here for thin F layers we observe a "double-kink" structure, see Fig. 7 (a) and(c). It corresponds to the DOS N f (E) with small minigap and finite subgap value smaller than unity. Such a DOS structure is typical in the presence of magnetic scattering and thin enough ferromagnetic interlayer, see Fig. 4 (a). The corresponding DOS energy dependencies are presented in the insets of the plots 6 (a

) [1/τ m ∆ = 3, h/∆ = 2, d f /ξ n = 0.5] and (b) [1/τ m ∆ = 3, h/∆ = 4, d f /ξ n = 0.5].
For finite magnetic scattering the non-monotonic features of CVC are smeared. We do not show the curves for d f /ξ n ≥ 2, since they do not significantly differ from the curves obtained from Eq. (29). temperatures T ≪ T C ). The value of this onset, according to Eq. ( 2), is ∆ + E g , where E g is the DOS minigap, 0 < E g < ∆.

Increasing exchange field, magnetic scattering and/or F layer thickness, the minigap closes and the current turns to zero at eV < ∆, having an onset at eV = ∆. The dependence of the minigap E g on the parameters characterizing the material properties of the ferromagnetic interlayer is discussed in Section III.

We conclude that we observe interesting features in the SIFS CVC if the DOS N f (E) near the insulating barrier has a nontrivial shape in the subgap region. In case when N f ≈ 1, these features disappear and the CVC coincide with those of the NIS tunnel junction, Eq. ( 28).

V. MACROSCOPIC QUANTUM TUNNELING IN A SIFS JUNCTION

In this section, motivated by experimental studies on the (thermal and quantum) switching 58 and quantum coherent oscillations 26 in SFS and SIFS junctions, we calculate the MQT escape rate in a current biased SIFS junction as shown in Fig. 8 (a). The CVC obtained in the previous section enable us to investigate the influence of the quasiparticle dissipation on MQT.

It is important to note that MQT can be used as a measurement process of a superconducting phase qubit. 59 Thus, the calculation of the MQT rate by taking into account the quasiparticle dissipation will be very important for analyzing the fidelity of measurement process for phase qubits. In the following calculation, for simplicity, we have ignored the influence from an environmental circuit on MQT which can be experimentally reduced by a noise filtering technique. 60 The partition function of a junction can be described by an imaginary-time functional integral over the macroscopic variable (the phase difference φ between two superconductors), [61][62][63][64] i.e.,

Z = Dφ (τ) exp - S eff [φ ] h . (30)
In the strong insulating barrier limit, i.e., γ B1 ≫ 1, the ef- fective action S eff is given by

S eff [φ ] = hβ 0 dτ M 2 ∂ φ (τ) ∂ τ 2 + U(φ ) + S α [φ ], (
α(τ) = h e ∞ 0 dω 2π exp(-ωτ) I V = hω e , (33)
at zero temperature. [61][62][63] As clearly seen from Figs. 6 and7, the CVC has a gap structure due to the isotropic superconducting gap in left superconductor electrode. In such a case, the dissipation kernel α decays exponentially as a function of imaginary time τ for |τ| ≫ h/∆. The typical dynamical time scale of the macroscopic variable φ is of the order of the inverse Josephson plasma frequency ω p = 2e|I c |/hC(1y 2 ) 1/4 which is much smaller than ∆. Thus, the phase varies slowly with the time scale given by h/∆, and we can expand φ (τ) -φ (τ ′ ) in Eq. (31b) about τ = τ ′ . This gives

S α [φ ] ≈ δC 2 hβ 0 dτ h 2e ∂ φ (τ) ∂ τ 2 , (34)
where

δC = 2 2e h 2 ∞ 0 dτα(τ)τ 2 . (35)
Hence, the dissipation action S α acts as a kinetic term so that the effect of the quasiparticles results in an increase of the capacitance, C → C + δC ≡ C ren .

In case of a thin ferromagnetic layer (d f = 0.5ξ n ) we numerically obtain δC ≈ h/∆R for CVC presented in Figs. 6 (a), 7 (a), (c) ["kink" and "double-kink" structures]. For thick ferromagnetic layer we can use Eq. ( 29) to calculate δC,

δC = 4h π∆R ∞ 0 dx x 2 ∞ 1 dz e -xz z 2 -1 ≈ 2h ∆R . (36)
For intermediate d f we numerically find δC ≈ (1 ÷ 2) h/∆R. Considering Nb as a superconductor (∆ = 1.3 meV) we therefore obtain δC ≈ (0.5 ÷ 1)r -1 pF, where r is the junction resistance R in Ω. To insure a small dissipative correction of capacitance, δC ≪ C we have a constraint, In order to see the effect of the quasiparticle dissipation on macroscopic quantum dynamics, we will investigate MQT in current-biased SIFS junctions. The MQT escape rate Γ from the metastable potential at zero temperature is given by 66

RC ≫ h/∆, (37)
Γ = lim β →∞ 2 β Im ln Z. ( 38
)
By using the Caldeira and Leggett theory, 68 the MQT rate is approximated as

Γ = ωp 2π √ 120πB exp(-B), (39)
where

ωp = 2e|I c | hC ren (1 -y 2 ) 1/4 , (40)
is the renormalized Josephson plasma frequency and B = S eff [φ B ]/h is the bounce exponent, that is the value of the action S eff evaluated along the bounce trajectory φ B (τ). The analytic expression for the bounce exponent is given by

B = 12 5e h 2e |I c |C ren 1 -y 2 5/4 . (41)
At high temperatures, the thermally activated decay dominates the escape process. Then the escape rate is given by the Kramers formula, Γ = ( ωp /2π) exp(-U 0 /k B T ), where U 0 is the barrier height. [START_REF] Weiss | Quantum Dissipative Systems[END_REF] Below the crossover temperature T * , the escape process is dominated by MQT. In the low dissipative (underdamping) cases, T * is approximately given by 66,67

T * = h ωp (y = y ) 2πk B , (42)
Here y = 1 0 dyP(y)y is the average switching current, where P(y) is the switching current distribution which is related to the escape rate Γ as 60

P(y) = 1 v Γ(y) exp - 1 v y 0 Γ(y ′ )dy ′ . (43)
In this equation, v ≡ |dy/dt| is the sweep rate of the external bias current. Importantly, T * is reduced in the presence of dissipative effects. 68 By using Eqs. ( 36) and ( 39), we calculate Γ and compare it with the case without the quasiparticle dissipation. In Fig. 9, we numerically plot Γ and Γ 0 for C = 800 pF and |I c | = 500 µA 13 and several values of R, where Γ 0 is the MQT escape rate in the absence of the quasiparticle dissipation [C ren → C in Eq. ( 39)]. As seen in this figure, Γ shows strong dependence on the junction resistance R, and Γ is almost identical to Γ 0 in the case of large R, e.g., R = 55 mΩ (≫ h/C∆ ≈ 0.63 mΩ) which corresponds to the actual SIFS junction. 13 We also calculate T * for a realistic case (R = 55 mΩ) and found that T * = 7.4mK for the dissipative case (C ren = C + δC) and T * = 7.5mK for the dissipation-less case (C ren = C). As expected, the T * suppression is small enough to allow experimental observations of MQT. Thus we can conclude that the influence of the quasiparticle dissipation on the macroscopic quantum dynamics of SIFS junction is very small for the case when the condition Eq. ( 37) is hold. This fact strongly suggests the great advantage of realistic SIFS junctions for qubit applications. The smallness of the quasiparticle dissipation in SIFS junctions is due to the superconducting gap in the left S electrode and the strong insulating barrier (γ B1 ≫ 1) between the left S and F layers.

It is important to note that such a weak quasiparticledissipation nature of MQT has been also predicted in π-junctions based on a S/ferromagnetic insulator (FI)/S junctions. 64,65 However no ferromagnetic insulator based Josephson junctions have been experimentally realized at present. On the other hands, the fabrication of SIFS junction is easily realized based on the current fabrication technology. [13][14][START_REF] Weides | Josephson Tunnel Junctions with Ferromagnetic Interlayer[END_REF]

VI. CONCLUSION

We have developed a quantitative theory, which describes the properties of the DOS and the current-voltage characteristics of a SIFS junction in the dirty limit. We considered the case of a strong insulating barrier in a SIFS junction such that the left S layer and the right F/S bilayer are decoupled. In this case we can obtain the current-voltage characteristics of a SIFS junction in the framework of standard tunnel theory. In order to calculate quasiparticle current we first calculated the DOS in the ferromagnetic layer of the F/S bilayer. We described the DOS behavior as a function of parameters characterizing properties of the ferromagnetic layer. In our theory we consider three such parameters: thickness of the ferromagnetic layer d f , exchange filed h, and magnetic scattering τ m . We have discussed the DOS properties paying special attention to the DOS oscillations at the Fermi energy. We have also proposed to measure the exchange field in experiments on weak ferromagnets by measuring the DOS peak at E = h. We compared the results, obtained with a self-consistent numerical method, with a known analytical DOS approximation, which is valid when the ferromagnetic layer is thick enough.

Using the numerically obtained DOS we have calculated the current-voltage characteristics of a SIFS junction and have observed features which are the signatures of the proximity effect in the S/F bilayer. We showed that there exists typical shape patterns of current-voltage characteristics related to the typical DOS structures in the ferromagnetic interlayer.

Finally, we have calculated the macroscopic quantum tunneling escape rate for the current biased SIFS junctions by taking into account the dissipative correction due to the quasiparticle tunneling. Based on this we concluded that the influence of the quasiparticle dissipation on the macroscopic quan-tum dynamics of SIFS junctions is small, which is a great advantage of SIFS junctions for qubit applications compared to other types of ferromagnetic π -junctions.

FIG. 2 :

 2 FIG. 2: (Color online) DOS N f (E) on the free boundary of the F layer in the F/S bilayer calculated numerically in the absence of magnetic scattering (1/τ m ∆ = 0) for different values of the F layer thickness d f , h/∆ = 4, T = 0.1T c . Parameters of the F/S interface are γ = γ B2 = 0.01. (a):d f /ξ n = 0.5, (b): d f /ξ n = 1, (c): d f /ξ n = 2, (d): d f /ξ n = 3.The approximate analytical solution, Eq. (15), is shown by dashed red lines.

FIG. 3 :

 3 FIG. 3: (Color online) DOS N f (E) on the free boundary of the F layer in the F/S bilayer calculated numerically in the absence of magnetic scattering (1/τ m ∆ = 0) for different values of the exchange field h. Parameters of the F/S interface are γ = γ B2 = 0.01, T = 0.1T c . Plots (a) and (b): d f /ξ n = 1; plots (c) and (d): d f /ξ n = 3. For plots (a) and (c) solid black line corresponds to h/∆ = 2, dashed red line to h/∆ = 2.5, dash-dotted blue line to h/∆ = 3. For plots (b) and (d) solid black line corresponds to h/∆ = 4, dashed red line to h/∆ = 5, dash-dotted blue line to h/∆ = 6.

FIG. 4 :

 4 FIG. 4: (Color online) DOS N f (E) at the free boundary of the F layer in the F/S bilayer calculated numerically for α m = 1/τ m ∆ = 0.5 (solid black line), α m = 1 (dashed red line), and α m = 3 (dash-dotted blue line) for different values of the F layer thickness d f , h = 4∆, T = 0.1T c . Parameters of the F/S interface are γ = γ B2 = 0.01. (a): d f /ξ n = 0.5, (b): d f /ξ n = 1, (c): d f /ξ n = 2, (d): d f /ξ n = 3. For plots (c) and (d) the curves with α m = 3 are not shown since they are of the order of unity at corresponding scale. Inset of the plot (a): N f (E) dependence for d f /ξ n = 0.5 for higher values of α m ; α m = 5 (solid black line), α m = 7 (dashed red line), α m = 15 (dash-dotted blue line).

FIG. 5 :

 5 FIG. 5: (Color online) Dependence of δ N f as a function of the F layer thickness d f for different exchange fields (a) and magnetic scattering times (b); dependence of δ N f as a function of exchange filed (c) and magnetic scattering time (d) for different d f . The temperature T = 0.1T C . Parameters of the F/S interface are γ = γ B2 = 0.01. (a): no magnetic scattering (1/τ m ∆ = 0), h/∆ = 2 (black solid line), h/∆ = 4 (red dashed line), h/∆ = 6 (blue dash-dotted line); (b): h/∆ = 4, α m = 1/τ m ∆ = 0 (black solid line), α m = 1 (red dashed line), α m = 3 (blue dash-dotted line); (c) no magnetic scat- tering, d f /ξ n = 0.5 (black solid line), d f /ξ n = 1 (red dashed line), d f /ξ n = 2 (blue dash-dotted line), d f /ξ n = 3 (black short-dashed line); (d) h/∆ = 4, d f /ξ n = 0.5 (black solid line), d f /ξ n = 1 (red dashed line), d f /ξ n = 2 (blue dash-dotted line), d f /ξ n = 3 (black short-dashed line).

FIG. 6 :

 6 FIG. 6: (Color online) Current-voltage characteristics of a SIFS junction in the absence of magnetic scattering for different values of the F-layer thickness d f . The temperature T = 0.1T c . The exchange field h = 0 (black line, which corresponds to the case of a SINS junction), h/∆ = 2 (blue dash-dotted line), and h/∆ = 4 (red short-dashed line). (a): d f /ξ n = 0.5; (b): d f /ξ n = 1, (c): d f /ξ n = 2, and (d): d f /ξ n = 3. Insets in (a), (b) and (c) are explained in the text.

FIG. 7 :

 7 FIG. 7: (Color online) Current-voltage characteristics of a SIFS junction for 1/τ m ∆ = 1 (black solid lines) and 1/τ m ∆ = 3 (red shortdashed lines) for different values of the F-layer thickness d f . The temperature T = 0.1T c . The exchange field h/∆ = 2: plots (a) and (b) and h/∆ = 4: plots (c) and (d). The thickness d f /ξ n = 0.5: plots (a) and (c) and d f /ξ n = 1: plots (b) and (d). Insets in (a) and (c) are explained in the text.

Figures 6 (

 6 a) and7 (a) show that for thin enough ferromagnetic layer the current has an onset in the interval [∆, 2∆] (for

FIG. 8 :

 8 FIG. 8: (Color online) (a) Schematic of a current biased SIFS Josephson junction. I ext is the external bias current. (b) Potential U(φ ) v.s. the phase difference φ between two superconductors. ω p is the Josephson plasma frequency of the junction.

  α(ττ ′ ) cos φ (τ) -φ (τ ′ ) 2 . (31b) In this equation, β = 1/k B T , M = C (h/2e) 2 is the mass (C is the capacitance of the junction) and the potential U(φ ) can be described by a tilted washboard potential [Fig. 8 (b)], i.e., U(φ ) = -E J [sgn(I c ) cos φ + yφ ] , (32) with y ≡ I ext /|I c |, where E J = h|I c |/2e is the Josephson coupling energy, I c is the Josephson critical current, and I ext is the external bias current, respectively. The dissipation kernel α(τ) is related to the quasiparticle current I under constant bias voltage V by

FIG. 9 :

 9 FIG. 9: (Color online) The MQT escape rate for a current-biased SIFS junction as a function of the bias current I ext for several values of the junction resistance R. Γ 0 (black solid line) and Γ (red dashed, yellow dot-dashed, and blue dotted lines) are the MQT escape rate without and with quasiparticle dissipation, respectively. Parameters are C = 800 pF, ∆ = 1.3 meV, and |I c | = 500 µA. 13i.e. the typical time constant RC of a SIFS junction should be much larger than the dynamical damping scale for the dissipation kernel α(τ). For example, in Ref.13 the following parameters of a Nb/Al 2 O 3 /Ni 0.6 Cu 0.4 /Nb SIFS junction were reported, C =800 pF and R =55 mΩ, which correspond to δC ≈ 10÷18 pF. Thus, even for the low resistive tunnel barrier in Ref.13, we have the condition Eq. (37) satisfied.In order to see the effect of the quasiparticle dissipation on macroscopic quantum dynamics, we will investigate MQT in current-biased SIFS junctions. The MQT escape rate Γ from the metastable potential at zero temperature is given by[START_REF] Weiss | Quantum Dissipative Systems[END_REF]

Acknowledgments

The authors thank A. I. Buzdin, E. Goldobin, D. A. Ivanov, A. V. Vedyayev, and M. Weides for useful discussions. This work was supported by NanoSCIERA project Nanofridge, ANR DYCOSMA, RFBR project N11-02-12065, JST-CREST, a Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan (Grant No. 22710096), the Invitation Program for Foreign Young Researchers in the G-COE program "Education and research center for material innovation", and the Spanish MICINN (Contract No. FIS2008-04209). F.S.B. thanks Intramural Special Project (Ref. 2009601036). A.A.G. thanks Dutch FOM for support. A.S.V. acknowledge the hospitality of Nanomaterials theory group, AIST, during his stay in Japan.