Bertrand Eynard
email: eynard@spht.saclay.cea.fr

Marcos Mariño
email: marcos@mail.cern.ch

Nicolas Orantin
email: orantin@spht.saclay.cea.fr

Holomorphic anomaly and matrix models

The genus g free energies of matrix models can be promoted to modular invariant, nonholomorphic amplitudes which only depend on the geometry of the classical spectral curve. We show that these non-holomorphic amplitudes satisfy the holomorphic anomaly equations of Bershadsky, Cecotti, Ooguri and Vafa. We derive as well holomorphic anomaly equations for the open string sector. These results provide evidence at all genera for the Dijkgraaf-Vafa conjecture relating matrix models to type B topological strings on certain local Calabi-Yau threefolds.

Introduction and conclusions

Topological string theory has been a fascinating laboratory to explore issues in string theory with important connections to other branches of physics and mathematics. The basic problem in topological string theory is to compute and understand closed and open string amplitudes on different geometric backgrounds. Particularly important among these are Calabi-Yau (CY) manifolds. Different techniques have been developed for the computation of these amplitudes. For type B topological strings on CY manifolds, a powerful method to solve the closed sector of the model are the holomorphic anomaly equations of [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF]. These equations control the t-dependence of the closed string amplitudes F (g) (t, t), and when combined with extra boundary conditions, they lead to explicit answers, see [START_REF] Huang | Topological string theory on compact Calabi-Yau: Modularity and boundary conditions[END_REF][START_REF] Grimm | Direct integration of the topological string[END_REF] for recent progress in this direction.

On the other hand, in some special backgrounds one can use large N dualities and geometric transitions to compute the holomorphic limit of F (g) (t, t), which will be denoted by F (g) (t). In a groundbreaking paper [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF], Dijkgraaf and Vafa conjectured that on certain noncompact Calabi-Yau manifolds, where the geometry reduces to a complex curve, the F (g) (t) are given by the genus g free energies of a matrix model in the 1/N expansion (see [START_REF] Mariño | Les Houches lectures on matrix models and topological strings[END_REF] for a review). The curve encoding the CY geometry is then identified as the classical spectral curve of the matrix model. Given the connection between topological strings and type II superstrings, this has made possible to compute certain protected quantities in supersymmetric gauge theories by using matrix model technology. The connection between topological strings on noncompact CY manifolds and matrix models has been extended recently to the mirrors of toric geometries [START_REF] Mariño | Open string amplitudes and large order behavior in topological string theory[END_REF].

The conjecture of Dijkgraaf and Vafa was verified at the planar level in [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF], and at genus one in [START_REF] Klemm | Gravitational corrections in supersymmetric gauge theory and matrix models[END_REF][START_REF] Dijkgraaf | Matrix models and gravitational corrections[END_REF]. In the simple case of the cubic matrix model, evidence for the conjecture was given at genus two in [START_REF] Huang | Holomorphic anomaly in gauge theories and matrix models[END_REF], where it was shown that the matrix model expression for F (g) (t) is the holomorphic limit of a particular solution to the holomorphic anomaly equations for the corresponding local curve.

In this paper we will give further evidence for the Dijkgraaf-Vafa conjecture and the related results of [START_REF] Mariño | Open string amplitudes and large order behavior in topological string theory[END_REF], by using recent advances in the solution of matrix models at all orders in the 1/N expansion obtained in [START_REF] Eynard | Topological expansion for the 1-hermitian matrix model correlation functions[END_REF][START_REF] Eynard | Topological expansion of the 2-matrix model correlation functions: Diagrammatic rules for a residue formula[END_REF][START_REF] Chekhov | Free energy topological expansion for the 2-matrix model[END_REF] and more recently in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. One of the outcomes of these advances is that, as explained in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], one can extrapolate the matrix model procedure and define a series of amplitudes F (g) (t) and correlation functions W g (p k ) for any algebraic curve H(x, y) = 0. When this curve is the classical spectral curve of a matrix model, one obtains in this way the 1/N expansion of the free energies and correlation functions. Using the results of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], we will show that it is possible to construct non-holomorphic free energies F (g) (t, t) and correlation functions which satisfy the following conditions:

• In the holomorphic limit t → ∞ they reduce to the holomorphic amplitudes F (g) (t) associated to the spectral curve.

• They are invariant with respect to the symplectic modular group of the curve. This is in contrast to their holomorphic limit, which does not have good modular properties [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

• They satisfy the holomorphic anomaly equations of [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF].

This gives a procedure to obtain the full non-holomorphic couplings F (g) (t, t) in the context of matrix models, by using the requirement of modular invariance as in [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF][START_REF] De Wit | Higher-order gravitational couplings and modular forms in N = 2, D = 4 heterotic string compactifications[END_REF], and shows that these couplings obey the equations of [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF] that characterize topological string amplitudes. This implies in particular that the matrix model free energy at genus g (extended non-holomorphically in this way) must be equal to the type B free energy of the noncompact Calabi-Yau manifolds considered in [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF][START_REF] Mariño | Open string amplitudes and large order behavior in topological string theory[END_REF], up to a holomorphic modular invariant quantity. Proving the Dijkgraaf-Vafa conjecture reduces now to proving that these holomorphic invariant functions vanish at all genera.

An interesting spinoff of our work is that we can determine the holomorphic anomaly equations for open string amplitudes with respect to the closed moduli t, at least in the local case. A direct string theory derivation of these equations has shown to be rather elusive [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF][START_REF] Antoniadis | Open string topological amplitudes and gaugino masses[END_REF], and we expect them to be useful in future investigations of the open sector in the general case. The equations we find turn out to be closely related to the holomorphic anomaly equation for correlation functions obtained in [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF], as we will explain in section 4.

This work has other interesting implications. If one considers the hyperelliptic curves that appear in Seiberg-Witten theory as classical spectral curves of a matrix model, the amplitudes F (g) constructed with the recursive procedure of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] should be closely related to the gravitational couplings of Seiberg-Witten theory introduced by Nekrasov in [START_REF] Nekrasov | Seiberg-Witten prepotential from instanton counting[END_REF]. It is known that these gravitational couplings can be promoted to non-holomorphic quantities that satisfy the holomorphic anomaly equations associated to the Seiberg-Witten curve [START_REF] Huang | Holomorphic anomaly in gauge theories and matrix models[END_REF][START_REF] Grimm | Direct integration of the topological string[END_REF]. It follows from the results in this paper that these gravitational couplings differ from the F (g) of the spectral curve at most in a holomorphic, modular invariant object, and it is natural to conjecture that they are in fact equal. The existence of a matrix model description of Nekrasov's F g has been proposed in [START_REF] Bertoldi | Large N gauge theories and topological cigars[END_REF].

Our work also leaves some open problems. First of all, although it goes a long way towards proving the conjecture in [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF], it would be nice to finish the proof by fixing the remaining ambiguity. At a more conceptual level, the t dependence of the amplitudes has been introduced, following [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], by requiring modular invariance and modifying the recursive construction of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] accordingly. It would be very interesting to find a more intrinsic way to introduce this non-holomorphic dependence in the theories which are defined by matrix integrals. This might shed some light on the problem of background independence in topological string theory and of the wavefunction behavior of the topological string partition function [START_REF] Witten | Quantum background independence in string theory[END_REF].

The organization of this paper is as follows. In section 2 we review the holomorphic anomaly equations and the connections between topological strings and matrix models. Most of the technology needed in this paper has been developed in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], and section 3 contains some results of that paper which will be essential in our derivation. Section 4 contains the derivation of the holomorphic anomaly equations. We actually provide a direct derivation and also a combinatorial derivation which makes contact with the recent results of [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF].

2 Holomorphic anomaly and topological strings

The holomorphic anomaly equations

Topological string theory on Calabi-Yau threefolds (see [START_REF] Mariño | Les Houches lectures on matrix models and topological strings[END_REF][START_REF] Neitzke | Topological strings and their physical applications[END_REF][START_REF] Vonk | A mini-course on topological strings[END_REF]) for recent reviews) is defined by coupling the twisted N = 2 sigma model to topological gravity and it comes in two versions: the A model (related to Gromov-Witten theory) and the B-model (related to deformations of complex structures). The topological amplitudes depend on a set of moduli that can be regarded as marginal deformations of the underlying N = 2 theory. If we denote the chiral and antichiral operators of the theory that correspond to marginal directions by φ I , φI , I = 1, • • • , n, we can deform the action in the standard way by using their descendants

S N =2 = S(0) + n I=1 t I Σg φ (2) I + n I=1 t I Σg φ (2) I (2.1)
In the closed string sector we are interested in computing the genus g free energy F (g) as well as the correlation functions of integrated chiral operators, C (g)

I 1 •••In = Σg φ (2) I 1 • • • Σg φ (2) In (2.2)
Since the t perturbation of the action can be written as a BRST-exact term with respect to the topological symmetry of the theory, one expects F (g) and the correlation functions to be t-independent. However, as discovered in [START_REF] Bershadsky | Holomorphic anomalies in topological field theories[END_REF][START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF], this is not the case once the twisted N = 2 sigma model is coupled to gravity. The BRST-exact terms give a boundary contribution that can be evaluated recursively, and one finds holomorphic anomaly equations for the t dependence of the amplitudes.

In order to write down the holomorphic anomaly equations, we need some basic ingredients from special geometry. The Zamolodchikov metric on the moduli space of marginal deformations is Kähler and can be derived from a Kähler function

G IJ = ∂ I ∂ J K.
(2.

3)

The amplitude at genus zero F 0 , also called the prepotential, leads to a holomorphic three point function

C IJK = ∂ 3 F 0 ∂t I ∂t J ∂t K .
(2.4)

The holomorphic two-point coupling is given by

τ IJ = ∂ 2 F 0 ∂t I ∂t J (2.5)
and plays an important role, as we will see later. Out of these quantities one can define a tensor

C IJ K = e 2K G IM G JN C M N K . (2.6)
The Zamolodchikov metric (2.3) defines a connection on moduli space with Christoffel symbols

Γ K IJ = G KM ∂ I G JM , Γ K IJ = G M K ∂ I G M J . (2.7)
It turns out that the F (g) are sections of a line bundle L 2-2g on the moduli space, and the covariant derivative acting on them is given by

D I = ∂ I -Γ I + (2 -2g)∂ I K, (2.8)
where K is the Kähler function appearing in (2.3). With these ingredients, the holomorphic anomaly equations of [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF] can be written as

∂ K F (g) = 1 2 C IJ K D I D J F g-1 + g-1 h=1 D I F h D J F g-h , g ≥ 2.
(2.9)

The equation for g = 1 is slightly different and reads

∂ I ∂ J F 1 = 1 2 C IKL C KL J - χ 24 -1 G IJ . (2.10)
Of course, the holomorphic anomaly equations determine F (g) (t, t) only up to a holomorphic function of the moduli. This integration constant is called the holomorphic ambiguity. The holomorphic limit of F (g) (t, t) is defined as

F (g) (t) = lim t→∞ F (g) (t, t). (2.11)
This is the limit which is for example appropriate to make contact with Gromov-Witten theory, in the A model. It has been observed [START_REF] De Wit | Higher-order gravitational couplings and modular forms in N = 2, D = 4 heterotic string compactifications[END_REF][START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF] that the non-holomorphic couplings F (g) (t, t) are invariant under the symmetry group acting naturally on the moduli space. In contrast, the holomorphic limit (2.11) breaks the modularity properties.

We finally point out that the correlation functions (2.2) can be written as covariant derivatives of the F (g) (t, t): C

(g)

I 1 •••In = D I 1 • • • D In F (g) , (2.12)
and they satisfy the holomorphic anomaly equation

∂ K C (g) I 1 •••In = 1 2 C M N K C (g-1) M N I 1 •••In + g r=0 n s=0 1 s!(n -s)! σ∈Sn C (r) M I σ(1) •••I σ(s) C (g-r) N I σ(s+1) •••I σ(n) -(2g -2 + n -1) n s=1 G KIs C (g) I 1 •••I s-1 I s+1 •••In .
(2.13)

The local Calabi-Yau case and the Dijkgraaf-Vafa conjecture

A particularly interesting class of topological strings are those defined on non-compact Calabi-Yau manifolds. These are commonly referred to as local Calabi-Yau manifolds.

In the local B model the Calabi-Yau geometry is encoded in a Riemann surface Σ n of genus n [START_REF] Klemm | Self-Dual Strings and N=2 Supersymmetric Field Theory[END_REF] (this should not be confused with the worldsheet of the string), and the resulting structure is significantly simpler than in the global case. The number of moduli is given by n, the genus of the curve, and they can be parametrized by

t I ∝ A I y(x)dx, I = 1, • • • , n, (2.14)
where A I are the compact A-cycles of Σ n , and y(x)dx is a meromorphic one-form on Σ n which can be constructed from the equation for the spectral curve. The proportionality constant in (2.14) will be determined below by a detailed comparison with the matrix model results. The imaginary part of the two-point coupling (2.5) is now a positive definite matrix, and the Zamolodchikov metric is simply given by

G I J = -i(τ -τ ) IJ . (2.15)
As noticed in [START_REF] Klemm | Local mirror symmetry at higher genus[END_REF][START_REF] Huang | Holomorphic anomaly in gauge theories and matrix models[END_REF][START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF], the holomorphic anomaly equations also simplify in the local case. First of all, the tensor C IJ K entering in (2.9) is now simply given by

C IJ K = G IM G JN C M N K = -[(τ -τ ) -1 ] IM [(τ -τ ) -1 ] JN C M N K . (2.16)
The covariant derivative acting on the F (g) amplitudes involves now only the Christoffel symbol (i.e. the term involving ∂ I K is no longer present). Notice that the Christoffel symbol can be written as

Γ L IJ = [(τ -τ ) -1 ] LM ∂ I τ JM , (2.17)
With these modifications, the holomorphic anomaly for F (g) , g ≥ 2 is still given by (2.9). For g = 1 the last term of (2.10) is not present. The holomorphic anomaly equations for the correlation functions are given by (2.13) but the contribution of the last line is absent in the local case. Since it will be useful later one, we now write the holomorphic anomaly equations for g ≥ 2, in the local case, in terms of a generating functional. We introduce,

Z = exp ∞ g=0 N 2-2g F (g) (t) , (2.18)
where N plays the role of a genus-counting parameter. We also define,

Z = exp -N 2 F (0) -F (1) Z, Ẑ = exp -N 2 F (0) Z. (2.19)
For g ≥ 2 we want to express the antiholomorphic derivative of Z in terms of Ẑ. Using that

∂ 2 IJ log Ẑ = 1 Ẑ ∂ 2 IJ Ẑ - 1 Z 2 ∂ I Ẑ∂ J Ẑ (2.20)
we can write (2.9) in the local case as

1 Z ∂ K Z = - 1 2N 2 1 Ẑ C IJ K ∂ 2 IJ Ẑ -Γ K IJ ∂ K Ẑ . (2.21)
In the local case, when the Calabi-Yau geometry reduces to a Riemann surface Σ n of genus n, the natural group acting on the moduli space is the mapping class group, and there is an induced action of the symplectic modular group

Sp(2n, Z) (2.22)
on the topological string amplitudes F (g) (t, t) and their holomorphic limit. The non-holomorphic amplitudes F (g) (t, t) of the local geometry turn out to be modular invariant under the action of this group [START_REF] Huang | Holomorphic anomaly in gauge theories and matrix models[END_REF][START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF][START_REF] De Wit | Higher-order gravitational couplings and modular forms in N = 2, D = 4 heterotic string compactifications[END_REF], while the F (g) (t) do not have good modular properties.

An interesting class of local Calabi-Yau backgrounds are the ones described by the equa-

tion uv = H(x, y), H(x, y) = y 2 -(W ′ (x)) 2 + f (x) (2.23)
where W (x), f (x) are polynomials of degree d + 1, d -1 respectively. The Riemann surface H(x, y) = 0 associated to this geometry is the hyperelliptic curve

y 2 = (W ′ (x)) 2 -f (x) (2.24)
of genus n = d -1. Dijkgraaf and Vafa conjectured in [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF] (see [START_REF] Mariño | Les Houches lectures on matrix models and topological strings[END_REF] for a review) that the holomorphic amplitudes F (g) (t) of type B topological string theory on the backgrounds of the form (2.23) are given by the 1/N expansion of a matrix model with spectral curve (2.24). This is in fact a matrix model with potential W (x) and generically with d -1 cuts around its extrema. This conjecture has been verified at the planar level in [START_REF] Dijkgraaf | Matrix models, topological strings, and supersymmetric gauge theories[END_REF], at genus one in [START_REF] Klemm | Gravitational corrections in supersymmetric gauge theory and matrix models[END_REF][START_REF] Dijkgraaf | Matrix models and gravitational corrections[END_REF], and at genus two in [START_REF] Huang | Holomorphic anomaly in gauge theories and matrix models[END_REF]. Another interesting class of local geometries is given by the mirrors of toric Calabi-Yau manifolds. In this case, the Riemann surface is a complex curve in "exponentiated" variables which can be written in the form

y 2 = M 2 (x) ℓ i=1 (x -x i ), (2.25)
where M(x) is in this case a transcendental function. As shown in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] and reviewed below, one can generalize the matrix model idea and associate a series of free energies F (g) (t) and correlation functions to any algebraic curve. It was recently proposed in [START_REF] Mariño | Open string amplitudes and large order behavior in topological string theory[END_REF] that holomorphic open and closed type B topological string amplitudes for the mirrors of toric background are given by the correlation functions and free energies, respectively, associated to the spectral curve (2.25) with the procedure of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

We then see that, for local CY backgrounds, where the geometry reduces to a complex curve, the holomorphic limit of topological string amplitudes in the B model has been conjectured to be given by the free energies and correlation functions associated to the complex curve by the recursive procedure of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. In the example (2.23) this coincides with the 1/N expansion of the free energy and correlators of a matrix model with potential W (x). In the rest of this paper we will give strong support to this conjecture by showing that the quantities defined in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] can be extended to non-holomorphic objects which have good modular properties and satisfy the holomorphic anomaly equation (2.9) in the local case.

3 Review of matrix models

Formal matrix models and algebraic geometry

Random matrix models became very popular in the 80's and 90's after it was discovered by 't Hooft [START_REF] Hooft | A planar diagram theory for strong interactions[END_REF] and then in [START_REF] Brézin | Planar Diagrams[END_REF] that the Feynman graphs of Hermitian matrix integrals are in fact discrete surfaces. Formal random matrix models are thus generating functions which count discrete maps of given topology. For instance the 1 matrix integral:

Z 1MM = dM e -N t Tr V (M ) (3.1)
counts discrete maps made of triangles, squares, ... k-gones where k ≤ deg V . The 2 matrix integral

Z 2MM = dM 1 dM 2 e -N t Tr (V 1 (M 1 )+V 2 (M 2 )-M 1 M 2 ) (3.2)
counts discrete maps made of polygons of two possible colors (or call it two possible spins ±). It was introduced by Kazakov [20] as an Ising model on a random map.

In all cases, the formal integral is defined as a formal series in its small t expansion, and it turns out that to any order in t, one has a "topological expansion":

F MM = ln Z MM = ∞ g=0 N 2-2g F (g) (3.3)
where F (g) is the generating function which counts maps of genus g. (see [START_REF] Eynard | Formal matrix integrals and combinatorics of maps[END_REF] for a review on this subject).

The goal was then to compute explicitly the coefficients F (g) in that expansion. Many methods have been invented, but the only one which has really succeeded so far is the so-called loop equation method which just amounts to integrate by parts or equivalently write Ward identities or Virasoro (or W-algebra) constraints. It was soon realized that loop equations allow in principle to find F (g) for any g [START_REF] Ambjorn | Matrix model calculations beyond the spherical limit[END_REF][START_REF] Akemann | Higher genus correlators for the Hermitian matrix model with multiple cuts[END_REF], but a more explicit solution (which can be represented diagrammatically) was found only recently first for the 1-matrix model [START_REF] Eynard | Topological expansion for the 1-hermitian matrix model correlation functions[END_REF], then for the 2-matrix model [START_REF] Eynard | Topological expansion of the 2-matrix model correlation functions: Diagrammatic rules for a residue formula[END_REF][START_REF] Chekhov | Free energy topological expansion for the 2-matrix model[END_REF], and more recently for the matrix model with external field [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

Then it was understood in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], that the solution of matrix models loop equation extends beyond the context of matrix models and is in fact a property of algebraic geometry in general. The diagrammatic method of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] allows to compute explicitly all the free energies F (g) appearing in (3.3) as well as all correlation functions (generating functions of discrete open surfaces of genus g with k boundaries):

Tr dx 1 x 1 -M . . . Tr dx k x k -M = ∞ g=0 N 2-2g-k W (g) k (x 1 , . . . x k ) (3.4)
where the mean value refers to the measure in (3.1) and the variables x i are fugacities associated to the lengths of boundaries.

Let us now give a brief summary of the construction of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

Given an arbitrary algebraic curve H(x, y) = 0, of genus n, one constructs recursively a sequence of multilinear symmetric meromorphic forms W (g) k and some scalars F (g) 's in terms of residues localized at branchpoints of the curve.

Many properties of those W (g) k 's and F (g) 's are described in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], in particular variations of conformal structure, integrable structure, invariance under symplectic transformations of the curve, homogeneity, and the one which interests us now is modular transformations.

Indeed, the construction of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] is based on the Bergmann kernel, which is defined for a given canonical basis of cycles. When one changes the basis of cycles, the Bergmann kernel changes, and the W (g) k 's and F (g) 's change accordingly. In [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], a "modified" Bergmann kernel was introduced in order to easily take into account those modular changes. The modified Bergmann kernel depends on an arbitrary complex symmetric matrix κ, and a modular change of cycles merely amounts (as far as only the Bergmann kernel is concerned) to a change of κ. Since the W (g) k 's and F (g) 's modular dependence is only in the Bergmann kernel, their modular properties are entirely encoded in their κ dependence. Thus in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], ∂κ and ∂F (g) ∂κ (3.5)

and they showed a striking similarity with holomorphic anomaly equations. The goal of this article is to prove that indeed one can derive the holomorphic anomaly equation from the construction of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF].

Variations of the matrix model free energies

Let us consider an algebraic curve H(x, y) of genus n and a canonical basis of cycles on it:

∀i, j = 1 . . . n , A i ∩ B j = δ ij , A i ∩ A j = 0 , B i ∩ B j = 0. (3.6)
There are n linearly independant holomorphic forms du i on H(x, y) normalized on the Acycles:

∀i, j = 1 . . . n ,

A i du j = δ ij (3.7)
and the Riemann matrix of period τ is a symmetric n × n matrix defined by

B j du i = τ ij . ( 3

.8)

Bergmann kernel

There exists a unique bilinear form B(p, q) with a unique double pole at p = q without residue and normalized on the A cycles:

B(p, q) ∼ p→q dz(p)dz(q) (z(p) -z(q)) 2 + finite , A B = 0. (3.9)
where z is any local parameter. B is often referred to as the Bergmann kernel in the literature. For instance, on a torus, B is the Weierstrass function.

Let us now introduce a new set of cycles, depending on an arbitrary complex symmetric matrix κ:

B := B -τ A , A := A -κB (3.10)
and define a κ-modified Bergmann kernel normalized on these new cycles by the constraints:

B(p, q) ∼ p→q dz(p)dz(q) (z(p) -z(q)) 2 + finite , A B = 0. (3.11)
This definition implies the relation:

B(p, q) = B(p, q) + 2iπ i,j du i (p) κ ij du j (q). (3.12)
Under a modular transformation:

τ → (C -τ D) -1 (τ E -F ) , E t C -F t D = Id , C t D = D t C , E t F = F t E (3.13)
the Bergmann kernel B changes, as well as the κ-modified Bergmann kernel, and the change is equivalent to a change of κ:

κ → D(C -τ D) -1 + (C -τ D) -1 t κ(C -τ D) -1 (3.14)
Therefore, the action of a modular transformation is equivalent to a change of κ, and all modular properties of the Bergmann kernel, and thus of the F (g) 's are encoded in the κ dependence.

One can verify that if κ = (ττ ) -1 , then the Bergmann kernel is modular invariant (see for instance a proof in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]).

Filling fractions

We also introduce the filling fractions

ǫ i := 1 2iπ A i ydx. (3.15)

Branch points

Let us also consider the branch points a i defined by dx(a i ) = 0. When a point p approaches a branch point a i , there exists a unique point p such that x(p) = x(p) and p → a i . Notice that p is defined locally near any branchpoint, but not globally (except for hyperelliptical curves of the form y 2 = P (x)).

The free energies and correlation functions are recursively defined as follows:

Correlation functions W (g) k = 0 if g < 0 (3.16) W (0) 1 (p) = 0 (3.17) W (0) 2 (p 1 , p 2 ) = B(p 1 , p 2 ) (3.18)
and

W (g) k+1 (p, p 1 , . . . , p k ) = i Res q→a i dE q (p) ω(q) g m=0 J⊂K W (m) j+1 (q, p J )W (g-m) k-j+1 (q, p K/J ) + W (g-1)
k+2 (q, q, p K )

(3.19)
with the notations (q being near a branchpoint):

ω(q) = (y(q) -y(q))dx(q), dE q (p) = 1 2 q q B(ξ, p) (3.20)
where the integration path lies entirely in a vicinity of a i . If J = {i 1 , i 2 , . . . , i j } is a set of indices, we write p J = {p i 1 , p i 2 , . . . , p i j }, and in the equation above we have K = {1, 2, . . . , k} and the summation over J is among all subsets of K.

It is important to notice that by construction, all W (g)

k 's have vanishing A-cycle integrals: A W (g) k = 0 (3.21)
which implies:

A W (g) k = κ B W (g) k (3.22)
Free energies.

For g > 1

F (g) = 1 2g -2 i Res q→a i Φ(q)W (g) 1 (q). (3.23)
where Φ(q) = q ydx is any antiderivative of ydx, i.e.

dΦ = ydx (3.24)

remark: it is proved in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] that

W (g) k = 1 2g + k -2 i Res a i ΦW (g)
k+1 .

(3.25) so that F (g) can be regarded as

F (g) = W (g) 0 .
Modular transformations and κ dependence.

From their definitions, it is clear that the only modular dependence, and also the only κ dependence of the W (g) k 's is through the Bergmann kernel, and therefore, a modular transformation is equivalent to a change of κ. Moreover, each W (g) k is a polynomial in κ of degree at most 3g -3 + 2k. In order to study modular transformations of the W (g) k 's, or their κ dependence, it is convenient to compute ∂W (g) k /∂κ. This was done in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] and the result is:

2iπ ∂ ∂κ ij W (g) k (p K ) = 1 2 r∈B j s∈B i W (g-1)
k+2 (p K , r, s)

+ 1 2 h L⊂K r∈B i W (h) |L|+1 (p L , r) s∈B j W (g-h) k-|L|+1 (p K/L , s) (3.26)
and in particular for k = 0:

2iπ ∂ ∂κ ij F (g) = 1 2 r∈B j s∈B i W (g-1) 2 (r, s) + 1 2 g-1 h=1 r∈B i W (h) 1 (r) s∈B j W (g-h) 1 (s) g ≥ 2.
(3.27) Also, notice that if κ = (τ -τ ) -1 , then all W (g) k 's and F (g) 's (with 2g +k ≥ 2) are modular invariant.

Variations wrt filling fractions

The variation of the Bergmann kernel wrt filling fractions follows from Rauch variational formula, and it was computed in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] that (in fact in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], κ was considered a constant, and thus we just have to add an extra ∂κ/∂ǫ term):

∂B(p, q) ∂ǫ = Res r→a B(r, p)B(r, q)du(r) dx(r)dy(r) -2iπ du t (p)κ ∂τ ∂ǫ κdu(q) + 2iπ du t (p) ∂κ ∂ǫ du(q) (3.28) which was written using an operator in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]:

∆ ǫ = ∂ ǫ + κ ∂τ ∂ǫ κ - ∂κ ∂ǫ ∂ ∂κ (3.29)
so that:

∆ ǫ B(p, q) = Res r→a B(r, p)B(r, q)du(r) dx(r)dy(r) (3.30)
This property was sufficient to ensure in [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF], that:

∆ ǫ W (g) k = - B W (g) k+1 (3.31)
4 Holomorphic anomaly equations and matrix models

In this section we give two proofs of the holomorphic anomaly equation for matrix models.

In the following, we consider the choice

κ = - 1 τ -τ (4.1)
in which case the F (g) 's and W

k 's are modular invariant [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] but not holomorphic.

A direct proof

We first consider the variation with respect to ǫ. The only antiholomorphic dependence of W (g) k on ǫ comes from its dependence on κ = 1/(ττ ), therefore

∂W (g) k ∂ǫ = ∂κ ∂ǫ ∂W (g) k ∂κ , ∂κ ∂ǫ = -κ ∂τ ∂ǫ κ , ∂ 2 F (0) ∂ǫ 2 = 2iπτ (4.2)
so that ∂W

(g) k ∂ǫ = - 1 2iπ κ ∂ 3 F (0) ∂ǫ 3 κ ∂W (g) k ∂κ . (4.3)
Then, using eq. (3.26):

∂W (g) k ∂ǫ = - 1 (2iπ) 2 κ ∂ 3 F (0) ∂ǫ 3 κ 1 2 r∈B j s∈B i W (g-1)
k+2 (p K , r, s)

+ h L⊂K W (h) |L|+1 (p L , r)W (g-h) k-|L|+1 (p K/L , s) (4.4)
Notice that with κ = 1/(ττ ), which satisfies:

κ ∂τ ∂ǫ κ = ∂κ ∂ǫ (4.5)
the differential operator ∆ of eq. (3.30) reduces to the usual derivative ∆ ǫ I = ∂/∂ǫ I , and therefore using eq. (3.31):

∂W (g) k ∂ǫ = - B W (g) k+1 = - B W (g) k+1 + τ A W (g) k+1 (4.6)
Acting again with ∂ ǫ we find (we use 3.22):

∂ 2 W (g) k ∂ǫ 2 = B B W (g) k+2 + ∂τ ∂ǫ A W (g) k+1 = B B W (g) k+2 + ∂τ ∂ǫ κ B W (g) k+1 (4.7)
and therefore:

∂ 2 W (g) k ∂ǫ 2 + ∂τ ∂ǫ κ ∂W (g) k ∂ǫ = B B W (g) k+2 (4.8)
Finally we get the equation:

∂W (g) k ∂ǫ = - 1 (2iπ) 2 κ ∂ 3 F (0) ∂ǫ 3 κ 1 2 ∂ 2 W (g-1) k ∂ǫ 2 + ∂τ ∂ǫ κ ∂W (g-1) k ∂ǫ + h L⊂K ∂W (h) l ∂ǫ ∂W (g-h) k-l ∂ǫ (4.9)
If we normalize the moduli as:

t I = (2iπ) 1 2 ǫ I . (4.10)
that last equation reads:

∂W (g) k ∂t = -κ ∂ 3 F (0) ∂t 3 κ 1 2 ∂ 2 W (g-1) k ∂t 2 + ∂τ ∂t κ ∂W (g-1) k ∂t + h L⊂K ∂W (h) l ∂t ∂W (g-h) k-l ∂t (4.11)
For the free energy, we set k = 0:

∂F (g) ∂t = -κ ∂ 3 F (0) ∂t 3 κ 1 2 ∂ 2 F (g-1) ∂t 2 + ∂τ ∂t κ ∂F (g-1) ∂t + g-1 h=1 ∂F (h) ∂t ∂F (g-h) ∂t (4.12) Since -κ IM ∂ 3 F (0) ∂t M ∂t N ∂t K κ JN = C IJ K , ∂τ JM ∂t I κ KM = -Γ K IJ (4.13)
we recover the holomorphic anomaly equation (2.9) for the free energies. The holomorphic anomaly equation for the correlation functions can be written (with the convention that

D I W (0) 1 = du I ) as ∂ K W (g) k = 1 2 C IJ K D I D J W (g-1) k + h L⊂K D I W (h) l D J W (g-h) k-l . (4.14)
In the context of topological string theory, these correlation functions correspond to open string amplitudes (see [START_REF] Mariño | Open string amplitudes and large order behavior in topological string theory[END_REF] for explicit results in the case of mirrors of toric CY manifolds), therefore the above equations can be regarded as holomorphic anomaly equations for open string amplitudes. These equations have been discussed in a general context in [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF][START_REF] Antoniadis | Open string topological amplitudes and gaugino masses[END_REF], and they are obtained by looking at the degenerations of Riemann surfaces with boundaries. The result (4.14) seems to indicates that the only degenerations that contribute in this case are the ones which involve closed string intermediate states. The first term in the r.h.s. of (4.14) corresponds to a node, while the sum in the second term corresponds to all the possible splittings of the Riemann surface Σ g,k into two Riemann surfaces Σ h,l and Σ g-h,k-l .

It is also easy to show that the correlation functions (2.12) are given by C (g)

I 1 •••I k = (-1) k (2iπ) k 2 B I 1 • • • B I k W (g) k (p K ) ( 4

A combinatorial proof

In the previous section we have given a simple proof of the holomorphic anomaly equations.

In [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF] (see also [START_REF] De Wit | Higher-order gravitational couplings and modular forms in N = 2, D = 4 heterotic string compactifications[END_REF]) it was noticed that these equations have a canonical solution of the form

F (g) (t, t) = F (g) (t) + Γ g ((Im τ ) -1 , ∂ I 1 • • • ∂ Im F r<g ), (4.18)
where Γ g is a polynomial functional in (Im τ ) -1 and in the derivatives of lower genera holomorphic free energies. In [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF], the coefficients of Γ g were interpreted diagrammatically as the possible degeneracies of a genus g Riemann surface. In this canonical solution, the antiholomorphic dependence of F (g) (t, t) enters only through Im τ . We will now show that this canonical solution can be obtained by iterating (3.26). This can be used to obtain a second proof of the holomorphic anomaly equations. Let F (g) (t, κ = 0) be the holomorphic free energy of the matrix model at genus g. When we turn on κ we get the κ-dependent function

d dκ g k = 1 2 1 2 + h g -h l k -l 2i -1 2i -1 2i 2i g -1 k
F (g) (t, κ) = 3g-3 m=0 κ m 1 m! d m dκ m F (g) (t, 0). ( 4.19)
The successive derivatives of F (g) (t, κ) at κ = 0 can be easily computed by iterating (3.26). We take into account that at κ = 0 we can substitute the contour integrals of the correlation functions W g k by ordinary (not covariant) derivatives of F (g) , as it follows from (4.7). It is convenient to represent (3.26) graphically, so that the terms obtained in the iteration are represented by graphs. This representation is shown in Fig. 1. Since the number of legs generated in this way is always even, we have labeled the connecting legs by 2i -1 and 2i. The result of the iteration is simply (this is easily seen from the graphical representation in Fig. 1):

F (g) (t, κ) = 3g-3 m=0 n I 1 ,...,I 2m =1 κ I 1 ,I 2 . . . κ I 2m-1 ,I 2m 1 2 m m! Gm A Gm , (4.20)
where G m is a connected degenerate surface with m propagators connecting 2m points labelled by 1, . . . , 2m, and in such a way that the point labeled by 2i -1 is connected by a propagator to the point labeled by 2i, i = 1, • • • , m. Each of these surfaces leads to a term A Gm . These terms are constructed as follows: if G m is made out of r Riemann surfaces of genera g 1 , • • • , g r with n 1 , • • • , n r punctures, respectively, (such that n i = 2m), and if the n i punctures of the i th Riemann surface are labeled J i = {J i,1 , . . . , J i,n i }, then

A Gm = r i=1 ∂ I J i,1 • • • ∂ I J i,n i F (g i ) (t, 0). (4.21)
an example of this procedure, we show in Fig. 2 the graphs that contribute to F 2 (t, κ).

Combinatorics

Notice because of the sum over indices I 1 , . . . , I 2m in (4.20), many graphs G m 's give the same contribution to F (g) (t, κ), and thus they better be computed only once with a multiplicity factor. The multiplicity factor is thus N Gm /(2 m m!), where N Gm is the number of inequivalent ways of relabeling the 2m punctures so that the puncture 2i -1 is linked by a propagator to 2i for all i = 1, . . . , m. Then notice that:

1 2 m m! = (2m -1)!! 2m! (4.22)
and notice that (2m -1)!! N Gm is the number of inequivalent ways of relabeling the 2m punctures (without any constraint on the labels of punctures which are linked together), i.e. it is the number of inequivalent pairings (compatible with G m ) of 2m points. This is exactly the multiplicity factor arising from Wick's theorem, i.e. F (g) (t, κ) is exactly the Feynman graph expansion of an integral whose propagator is N -2 κ and whose vertices are the N

2-2g i ∂ I J i,1 • • • ∂ I J i,n i F (g i ) (t, 0).
Notice that because we have chosen W

= 0 and B B W (0) 2 = 0, A m vanishes if G m contains surfaces of genus zero with less than 3 punctures. Therefore there is no vertices corresponding to the first and second derivative of F (0) .

Notice also that the recursive application of Fig. 1 generates only connected graphs G m 's, and thus we are in fact computing the Feynman graph expansion of the log of an integral. Therefore

F (t, κ) = ∞ g=0
N 2-2g F (g) (t, κ) = ln (Z(t, κ)) (4.23) Figure 2: The graphs that contribute to F 2 (t, κ) after iterating Fig. 1.

+ + + + + + κ 2 2 2 • 2!
is exactly the Feynman graph expansion the following integral

Z(t, κ) = dη exp - 1 2 N 2 (η -t)κ -1 (η -t) + F (η, 0) -N 2 (η -t) I ∂ I F (0) (t, 0) - 1 2 N 2 (η -t) I ∂ 2 IJ F (0) (t, 0)(η -t) J . (4.24)
Indeed, the Feynman graph expansion, is nothing but the saddle point approximation, and here the saddle point equation to large N leading order is κ -1 IJ (ηt) J -∂ I F (0) (η, 0) + ∂ I F (0) (t, 0) + ∂ 2 IJ F (0) (t, 0)(ηt) J = 0, which is clearly satisfied at η = t. The propagator is simply κ IJ N 2 (4.25)

and the vertices are given by

N 2-2g k! ∂ I 1 • • • ∂ I k F (g) (t), (4.26)
where k ≥ 3 for g = 0 and k ≥ 1 for g ≥ 1.

We will now show that, for the choice κ = (ττ ) -1 , the generating functional (4.24) satisfies the holomorphic anomaly equations of [START_REF] Bershadsky | Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes[END_REF] in the local case. We first note that there is a cancellation, Z(t, (ττ ) -1 ) = dη exp -1 2 N 2 (ηt)τ (ηt) + F (η, 0) -N 2 (ηt) I ∂ I F (0) (t, 0) (4.27)

We compute now:

1 Z ∂ I Z = N 2 (τ IJ -τ IJ ) (η -t) J + N 2 ∂ I F (0) , (4.28)
where denotes an average in the integral (4.27). It follows that, On the other hand, we compute: we finally find that

1 Z ∂ K Z = -
N 2 Z ∂ K Z = 1 2 Ẑ C IJ K ∂ 2 IJ Ẑ -Γ L IJ ∂ L Ẑ , (4.36)
where the tensor C IJ K and the Christoffel symbol Γ L IJ are given by (2.16) and (2.17), respectively. In other words we obtain the holomorphic anomaly equations (2.21) which is equivalent to (2.9).

The expansion in [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF] is precisely the one that is obtained by expanding the integral (4.24). This means that the iteration of Fig. 1 leads to a representation of F (g) (t, κ) which agrees with the canonical solution of [START_REF] Aganagic | Topological strings and (almost) modular forms[END_REF] for the choice (4.1).

D

  .[START_REF] Eynard | Topological expansion of the 2-matrix model correlation functions: Diagrammatic rules for a residue formula[END_REF] since the covariant derivative D I = ∂ I -Γ I commutes with B when acting on a differential f (p)dp with vanishing A-cycle integrals:I f (p)dp -Γ I B f (p)dp -∂ I τ A f (p)dp = B D I f (p)dp -Γ I B f (p)dp -∂ I τ κ B f (p)dp = B D I f (p)dp (4.16)Note that the last equality requires explicitly the conditionA f (p)dp = A f (p)dpκ B f (p)dp = 0. (4.17) If we now integrate (4.14) over the cycles B I 1 , • • • , B I k , one finds the holomorphic anomaly equation for the correlation function (2.13) in the local case (i.e. without the last line). This indicates again that the open string anomaly equations (4.14) involve only the closed string sector in the intermediate states.

Figure 1 :

 1 Figure 1: A graphic representation of the equation (3.26).

1 Ẑ 2

 12 ∂ I Ẑ = N 2 (τ IJτ IJ ) (ηt) J , JI Ẑ = -N 2 (τ IJ -τ IJ )-N 2 ∂ J τ IK (η-t) K +N 4 (τ -τ ) IK (τ -τ ) JL (η-t) K (η-t) L .(4.30)

1 2 N 2 CF ( 1 )

 1221 IJK (ηt) I (ηt) that the non-holomorphic F(1) (κ) generated by the integral (4.27) is (t, κ) = F(1) (t) + ∆(κ). (4.34)We express now the quadratic average in(4.31) in term of the double derivative (4.30). If we use ∂ K ∆ = κ IJ ∂ K τIJ , (4.35)

Acknowledgments

We would like to thank Mina Aganagic, Robert Dijkgraaf, Amir Kashani-Poor, Albrecht Klemm, Andrei Okounkov and Pierre Vanhove for useful and fruitful discussions on this subject. This work is partly supported by the Enigma European network MRT-CT-2004-5652, by the ANR project Géométrie et intégrabilité en physique mathématique ANR-05-BLAN-0029-01, by the Enrage European network MRTN-CT-2004-005616, by the European Science Foundation through the Misgam program, by the French and Japaneese governments through PAI Sakurav, by the Quebec government with the FQRNT.

the following quantities were computed:


Marine Schott

Anne Bernand-Mantel
email: anne.bernand-mantel@neel.cnrs.fr

Laurent Ranno

Stéfania Pizzini

Jan Vogel

Hélène Béa

Claire Baraduc

Stéphane Auffret

Gilles Gaudin

Dominique Givord

The skyrmion switch: turning magnetic skyrmion bubbles on and off with an electric field

Keywords: room temperature, magnetic switching, skyrmions, electric field, thin films, spintronics, perpendicular magnetic anisotropy

Nanoscale magnetic skyrmions are considered as potential information carriers for future spintronics memory and logic devices. Such applications will require the control of their local creation and annihilation, which involves so far solutions that are either energy consuming or difficult to integrate. Here we demonstrate the control of skyrmion bubbles nucleation and annihilation using electric field gating, an easily integrable and potentially energetically efficient solution. We present a detailed stability diagram of the skyrmion bubbles in a Pt/Co/oxide trilayer and show that their stability can be controlled via an applied electric field. An analytical bubble model, with the Dzyaloshinskii-Moriya interaction imbedded in the domain wall energy, account for the observed electrical skyrmion switching effect. This allows us to unveil the origin of the electrical control of skyrmions stability and to show that both magnetic dipolar interaction and the Dzyaloshinskii-Moriya interaction play an important role in the skyrmion bubble stabilization.

stabilized by the magnetic dipolar energy, a skyrmion bubble is stabilized by both dipolar and Dzyaloshinskii-Moriya interaction (DMI) energies. A second difference between classical bubbles and skyrmionic bubbles is their homo-chirality which allows their current driven motion via spin-orbit torques 7,9,13,[START_REF] Yu | Room-Temperature Creation and Spin-Orbit Torque Manipulation of Skyrmions in Thin Films with Engineered Asymmetry[END_REF] . This discovery has triggered a growing interest for the use of magnetic skyrmions as building blocks for memories [START_REF] Kiselev | Chiral Skyrmions in Thin Magnetic Films: New Objects for Magnetic Storage Technologies[END_REF][START_REF] Fert | Skyrmions on the Track[END_REF][START_REF] Sampaio | Stability and Current-Induced Motion of Isolated Magnetic Skyrmions in Nanostructures[END_REF][START_REF] Zhang | Skyrmion-Skyrmion and Skyrmion-Edge Repulsions in Skyrmion-Based Racetrack Memory[END_REF][START_REF] Zhang | Magnetic Skyrmion Logic Gates: Conversion, Duplication and Merging of Skyrmions[END_REF][START_REF] Kang | Voltage Controlled Magnetic Skyrmion Motion for Racetrack Memory[END_REF][START_REF] Wiesendanger | Nanoscale Magnetic Skyrmions in Metallic Films and Multilayers: A New Twist for Spintronics[END_REF] which is reminiscent of the extensive research on magnetic bubble memories in the 1970s [START_REF] Bobeck | Properties and Device Applications of Magnetic Domains in Orthoferrites[END_REF] . Such applications will require creating, deleting and moving skyrmions. The recent observation of current induced displacement of skyrmion bubbles using relatively low current densities reported at RT 9,13,[START_REF] Yu | Room-Temperature Creation and Spin-Orbit Torque Manipulation of Skyrmions in Thin Films with Engineered Asymmetry[END_REF][START_REF] Yu | Skyrmion Flow near Room Temperature in an Ultralow Current Density[END_REF] is very promising. The creation and annihilation of skyrmions have been addressed theoretically [START_REF] Sampaio | Stability and Current-Induced Motion of Isolated Magnetic Skyrmions in Nanostructures[END_REF][START_REF] Koshibae | Creation of Skyrmions and Antiskyrmions by Local Heating[END_REF][START_REF] Zhou | Dynamically Stabilized Magnetic Skyrmions[END_REF] and experimentally 6,[START_REF] Hsu | Electric-Field-Driven Switching of Individual Magnetic Skyrmions[END_REF][START_REF] Nii | Uniaxial Stress Control of Skyrmion Phase[END_REF][START_REF] Finazzi | Laser-Induced Magnetic Nanostructures with Tunable Topological Properties[END_REF] by different techniques such as spin transfer torque 6,[START_REF] Sampaio | Stability and Current-Induced Motion of Isolated Magnetic Skyrmions in Nanostructures[END_REF][START_REF] Zhou | Dynamically Stabilized Magnetic Skyrmions[END_REF] , heat [START_REF] Koshibae | Creation of Skyrmions and Antiskyrmions by Local Heating[END_REF][START_REF] Finazzi | Laser-Induced Magnetic Nanostructures with Tunable Topological Properties[END_REF] and strain [START_REF] Nii | Uniaxial Stress Control of Skyrmion Phase[END_REF] . However, those techniques are either energy consuming or difficult to integrate in functional devices. In analogy to what has been developed to switch small ferromagnetic elements [START_REF] Wang | Voltage-Induced Switching in Magnetic Tunnel Junctions with Perpendicular Magnetic Anisotropy[END_REF] , the use of electric field gating to manipulate skyrmions offers several advantages: the low power consumption, the possibility to act locally and an easy integration. Furthermore, a major issue of current induced nucleation, namely, unwanted current displacement of skyrmions during the writing process, is naturally absent using electric field writing. Those advantages have motivated some theoretical works concerning the electric control of skyrmions [START_REF] Kang | Voltage Controlled Magnetic Skyrmion Motion for Racetrack Memory[END_REF][START_REF] Fook | Gateable Skyrmion Transport via Field-Induced Potential Barrier Modulation[END_REF][START_REF] Upadhyaya | Electric-Field Guiding of Magnetic Skyrmions[END_REF][START_REF] Nakatani | Electric Field Control of Skyrmions in Magnetic Nanodisks[END_REF] . The only experimental attempt for skyrmion electrical switching has been carried out in an epitaxial material and at low temperature [START_REF] Hsu | Electric-Field-Driven Switching of Individual Magnetic Skyrmions[END_REF] . In this letter, we report on the observation of thermally activated nucleation and annihilation of magnetic skyrmion bubbles at RT and we demonstrate the possibility to electrically switch them on and off under a constant biasing magnetic field.

We have selected for this study a 𝑃𝑡/𝐶𝑜/𝑜𝑥𝑖𝑑𝑒 trilayer system where large DMI energy values have been reported [START_REF] Pizzini | Chirality-Induced Asymmetric Magnetic Nucleation in Pt/Co/AlOx Ultrathin Microstructures[END_REF][START_REF] Vaňatka | Velocity Asymmetry of Dzyaloshinskii Domain Walls in the Creep and Flow Regimes[END_REF][START_REF] Tetienne | The Nature of Domain Walls in Ultrathin Ferromagnets Revealed by Scanning Nanomagnetometry[END_REF][START_REF] Belmeguenai | Interfacial Dzyaloshinskii-Moriya Interaction in Perpendicularly Magnetized Pt/Co/AlO Ultrathin Films Measured by Brillouin Light Spectroscopy[END_REF] . In addition, we have chosen thicknesses for which the 𝐶𝑜 layer is close to the ferromagnetic-paramagnetic transition at RT, as enhanced electric field effects have been observed in this conditions [START_REF] Chiba | Electrical Control of the Ferromagnetic Phase Transition in Cobalt at Room Temperature[END_REF] . The 𝐶𝑜 layer presents a slight thickness gradient (2%/mm) created by the oxidation of a wedge shaped 𝐴𝑙 top layer (see Figure 1b). A high-k 𝐻𝑓𝑂 2 dielectric layer and a transparent top Indium Tin oxide (𝐼𝑇𝑂) electrode have been deposited on top of the 𝑃𝑡/𝐶𝑜/𝑜𝑥𝑖𝑑𝑒 trilayer, allowing magnetic characterization using magneto-optical Kerr effect under applied electric fields.

The magnetic domains are imaged at RT using a Kerr microscope in polar geometry. At low magnetic field a spontaneously demagnetized labyrinthine state is observed (Figure 1a). The formation of labyrinthine domains is the consequence of a balance between competing energies: the magnetic dipolar energy and the domain wall (DW) energy. A variation of the labyrinthine domain width 𝐿 is observed along the wedge direction (Figure 1c). To understand this variation we have used an analytical model [START_REF] Hubert | Magnetic Domains, The Analysis of Magnetic Microstructures[END_REF][START_REF] Kaplan | The Domain Structure in Ultrathin Magnetic Films[END_REF] describing for ultrathin films the labyrinthine domain width 𝐿 = 𝛼. 𝑡. exp(𝜋𝐿 0 /𝑡) as function of the magnetic film thickness 𝑡, the characteristic dipolar length 𝐿 0 = 𝜎 𝑊 /𝜇 0 . 𝑀 𝑠 Using the experimental effective anisotropy value 𝐾 𝑒𝑓𝑓 = 0.30 MJ/m 3 (see supplementary S4) and an exchange constant 𝐴 𝑒𝑥 in the 5 -15 pJ/m range we find 𝜎 𝑤 𝐵𝑙𝑜𝑐ℎ = 5 -8 mJ/m² which is much larger than the experimental DW energy. The deviation from the Bloch formula can be explained by the presence of interfacial DMI in the trilayer which reduces 𝜎 𝑊 according to the expression 13,[START_REF] Heide | Dzyaloshinskii-Moriya Interaction Accounting for the Orientation of Magnetic Domains in Ultrathin Films: Fe/W(110)[END_REF] : 𝜎 𝑊 = 4 √ 𝐴 𝑒𝑥 𝐾 𝑒𝑓𝑓 -𝜋𝐷 where D is the DMI energy in mJ/m². We deduce a DMI term 𝐷 = 1 -2 mJ/m² which is consistent with the values obtained in similar 𝑃𝑡/𝐶𝑜/𝑜𝑥𝑖𝑑𝑒 trilayers [START_REF] Pizzini | Chirality-Induced Asymmetric Magnetic Nucleation in Pt/Co/AlOx Ultrathin Microstructures[END_REF][START_REF] Vaňatka | Velocity Asymmetry of Dzyaloshinskii Domain Walls in the Creep and Flow Regimes[END_REF][START_REF] Tetienne | The Nature of Domain Walls in Ultrathin Ferromagnets Revealed by Scanning Nanomagnetometry[END_REF][START_REF] Belmeguenai | Interfacial Dzyaloshinskii-Moriya Interaction in Perpendicularly Magnetized Pt/Co/AlO Ultrathin Films Measured by Brillouin Light Spectroscopy[END_REF] . The low DW energy in our system is at the origin of the observed spontaneous demagnetization of the layer at RT on a few second timescale. This thermally activated demagnetization is due to a combination of spontaneous nucleation of magnetic domains and thermally activated motion of domain walls under zero applied magnetic field (see supplementary S5). Such thermal movement of labyrinthine magnetic domains has been reported

previously in ultrathin films where the DW energy was lowered by the proximity of the spin reorientation transition [START_REF] Won | Magnetic Stripe Melting at the Spin Reorientation Transition in Fe Ni Cu(001)[END_REF][START_REF] Bergeard | Dynamic Fluctuations and Two-Dimensional Melting at the Spin Reorientation Transition[END_REF] . This DW mobility, in combination with the high nucleation rate, allows reaching an equilibrium demagnetized state in a few seconds in the major part of the image of Figure 1a (between 𝑝 0 and 𝑝 3 ). Above 𝑝 0 , the demagnetized state starts to deviate from the equilibrium one and we decided not to extract the domain width in this region. We now discuss the origin of the decrease by a factor of 3 of the domain width 𝐿 along the wedge (Figure 1c). The 𝐶𝑜 thickness (𝑡) variation is of the order of 1% along the considered region and its effect on 𝐿 in the expression 𝐿 = 𝛼. 𝑡. 𝑒𝑥𝑝( 𝜋𝜎 𝑊 𝜇 0 .𝑀 𝑠 2 /𝑡) is negligible.

However, as 𝑡 decreases, a reduction of 𝑀 𝑠 by 8.4% is observed between 𝑝 0 and 𝑝 3 due to a variation of 𝑇 𝐶 with 𝑡 and the proximity of 𝑇 𝐶 at the observation temperature [START_REF] Parakkat | Tailoring Curie Temperature and Magnetic Anisotropy in Ultrathin Pt/Co/Pt Films[END_REF] (see supplementary S2). According to its analytical expression, the domain width 𝐿 should increase when 𝑀 𝑠 , and thus the magnetic dipolar energy, decreases. However, the opposite is observed in Figure 1c. We deduce that a strong decrease of the DW energy 𝜎 𝑊 , over-compensating the 𝑀 𝑆 variation is at the origin of the decrease of 𝐿. If we assume a constant DMI energy term, the variation of the DW energy 𝜎 𝑊 = 4 √ 𝐴 𝑒𝑥 𝐾 𝑒𝑓𝑓 -𝜋𝐷 must be related to 𝐴 𝑒𝑥 and 𝐾 𝑒𝑓𝑓 variations. This is confirmed by our experimental observations: a reduction of 𝐾 𝑒𝑓𝑓 with 𝑡 is measured (see supplementary S3) and a diminution of 𝐴 𝑒𝑥 is expected from the 𝑀 𝑆 decrease with 𝑡 [START_REF] Eyrich | Effects of Substitution on the Exchange Stiffness and Magnetization of Co Films[END_REF][START_REF] Moreno | Temperature-Dependent Exchange Stiffness and Domain Wall Width in Co[END_REF] .

When a perpendicular magnetic field larger than 0.1 mT is applied, micron-sized bubbles appear (Figure 1a). This transition from a labyrinthine pattern into a bubble lattice with applied magnetic field is unusual for a classical bubble system. In the classical case the magnetic domains would shrink in width and length but the number of magnetic objects would remain constant [START_REF] Hubert | Magnetic Domains, The Analysis of Magnetic Microstructures[END_REF] , leaving very few bubbles. This situation is observed on the right side of the image in Figure 1a. As we go from right to left, the increasing density of bubbles domains indicates the presence of thermally activated nucleation. This thermally activated process is promoted by the reduction of the DW energy 𝜎 𝑤 , due to the presence of DMI. This implies a skyrmionic nature of the nucleated bubbles, which is confirmed by the observation of their unidirectional motion against electron flow, as expected for this system (see supplementary S6). As the nucleation of a skyrmion bubble by coherent rotation of a 1 µm domain is too high in energy (𝐾𝑉 > 10 4 𝑘 𝐵 𝑇 293𝐾 ) the mechanism for the observed skyrmion bubbles creation has to involve expansion from a nanoscale skyrmion. We can estimate the nucleation energy of a skyrmion bubble when its diameter is reduced to the DW width 𝛿 𝑤 by 𝜎 𝑤 . 𝑡. 2. 𝜋. 𝛿 𝑤 ~10𝑘 𝐵 𝑇 293𝐾 . Following a Néel Brown model, the RT nucleation rate of such skyrmion bubbles is high (~0.5 MHz for 𝜏 0 = 0.1 𝑛s). To estimate whether the nucleated skyrmion bubble is expected to annihilate or to grow and reach its equilibrium diameter, we have calculated analytically the energy of an isolated skyrmion bubble using a "thin wall" model which is valid for negligible DW width (𝛿 𝑤 ~10 nm in our case) compared to the bubble diameter (few 𝜇𝑚) (see Methods). The skyrmion bubble energy 𝐸 𝑠𝑏 (𝑡, 𝑀 𝑠 , 𝜎 𝑤 ) is the sum of DW, Zeeman and dipolar terms. The energy of an individual bubble relative to the saturated state can be written as: where 𝜎 𝑤 is the DW energy, 𝑡 the layer thickness, 𝑅 the bubble radius, 𝑀 𝑠 the saturation magnetization, 𝜇 0 𝐻 the applied magnetic field, and 𝐼(𝑑) the stray field energy gain where 𝑑 = 2𝑅/𝑡 (see Methods). The skyrmion bubble energy ∆𝐸 𝑠𝑏 versus the bubble diameter (Figure 2.g), obtained using the experimental 𝑡, 𝑀 𝑠 and 𝜎 𝑤 parameters at positon 𝑝 1 , presents two local extrema. The existence of a local energy minimum in the model is consistent with the experimental observation of an equilibrium diameter for the skyrmion bubbles. To reach this equilibrium diameter, a nucleated nanoscale skyrmion will have to grow and overcome a nucleation energy barrier 𝐸 𝑛 , which corresponds to the local maximum. We define the annihilation energy barrier 𝐸 𝑎 as the difference between the local maximum and minimum (see Figure 2.g), assuming that the lower energy path for annihilation is the compression [START_REF] Bessarab | Method for Finding Mechanism and Activation Energy of Magnetic Transitions, Applied to Skyrmion and Antivortex Annihilation[END_REF][START_REF] Rohart | Path to Collapse for an Isolated Neel Skyrmion[END_REF][START_REF] Lobanov | Mechanism and Activation Energy of Magnetic Skyrmion Annihilation Obtained from Minimum Energy Path Calculations[END_REF][START_REF] Cortés-Ortuño | Thermal Stability and Topological Protection of Skyrmions in Nanotracks[END_REF] sites is observed at fixed electric field. The variation of skyrmion bubbles density with the electric field has been characterized by image processing (see Figure 2f). To explain the observed behavior, the analytical skyrmion bubble model described previously is used. The 𝑡, 𝑀 𝑠 and 𝜎 𝑤 parameters and their variations with the electric field, used as an input for the simulation, are determined experimentally (see supplementary S3 and S8). When the voltage is switched from -5 V to +5 V, 𝑀 𝑠 and 𝜎 𝑤 vary by respectively +0.1 MA/m and +0.32 mJ/m². These experimental parameters are used in the simulation in Figure 2g where the energy difference between a ferromagnetic state and an isolated bubble is plotted as a function of the bubble diameter. The observed decrease in the skyrmion bubbles number with increasing electric field is consistent with the increase in the nucleation energy barrier 𝐸 𝑛 combined with a decrease of the annihilation energy barrier 𝐸 𝑎 obtained in the simulation. In Figure 2 d,e differential images are shown where the blue/red bubbles correspond to bubbles which have appeared/disappeared. In Figure 2d, both blue and red bubbles are observed. This is explained by a domination of nucleation processes as 𝐸 𝑛 ≤ 23 𝑘 𝐵 𝑇 293𝐾 and 𝐸 𝑎 ≥ 23 𝑘 𝐵 𝑇 293𝐾 in the -5 V to 0 V range (Figure 2g). The observation of both blue and ref bubbles in that voltage range is a consequence of the displacement of existing bubbles due to a continuous nucleation of bubbles. On the contrary the 0 𝑉 to +5 𝑉 range is dominated by annihilation processes with 𝐸 𝑎 ≤ 23 𝑘 𝐵 𝑇 293𝐾 and 𝐸 𝑛 ≥ 23 𝑘 𝐵 𝑇 293𝐾 as confirmed by a large majority of red bubbles in Figure 2e.

∆𝐸 𝑠𝑏 (𝑅,
Now that the electric field control of skyrmion bubble density is demonstrated, we carry out the proof of concept of an electrical skyrmion switch by turning the electric field to higher voltages. In Figure 3a-b, we show the electric field switching from a state with a very high skyrmion bubbles density to a state with no skyrmion bubbles. This electric switching of skyrmion bubbles is reversible and reproducible sequentially as we can see in Figure 3d where the average Kerr intensity is recorded below the electrode as function of time while the voltage is switched between +/-20 V (see supplementary S9). To understand the origin of this switching effect, polar Kerr hysteresis loops with in plane applied magnetic field have been measured to estimate the variation of 𝑀 𝑆 and 𝐾 𝑒𝑓𝑓 under electric field (see supplementary S8). When the voltage varies from +20 V to -20 V, 𝑀 𝑆 and 𝐾 𝑒𝑓𝑓 are modified respectively by -47% and -64%. We deduce from 𝐾 𝑒𝑓𝑓 = 𝐾 𝑠 /𝑡 -1/2𝜇 0 𝑀 𝑆 2 a 𝐾 𝑠 variation ∆𝐾 𝑠 /∆𝐸 = 830 fJ/V/m. The electric field variations of the magnetic properties is most likely due to a modification of the electron density of state in the 𝐶𝑜 and enhanced by the proximity of 𝑇 𝐶 as observed in previous works on 𝑃𝑡/𝐶𝑜/𝑜𝑥𝑖𝑑𝑒 in the ultrathin regime [START_REF] Chiba | Electrical Control of the Ferromagnetic Phase Transition in Cobalt at Room Temperature[END_REF] . A variation of 𝐷 with electric field may also occur in our experiment as predicted theoretically [START_REF] Yang | Controlling Dzyaloshinskii-Moriya Interaction via Chirality Dependent Layer Stacking, Insulator Capping and Electric Field[END_REF] . It is possible to extract D from the expression 𝜎 𝑊 = 4 √ 𝐴 𝑒𝑥 𝐾 𝑒𝑓𝑓 -𝜋𝐷 using the experimental values of 𝐾 𝑒𝑓𝑓 and 𝜎 𝑊 . However, as 𝑀 𝑠 is varying with the electric field (due to the variation of 𝑇 𝐶 ) the exchange constant 𝐴 𝑒𝑥 is also expected to vary and the detailed analysis necessary to separate the D variation from A ex is beyond the scope of this study. The variation of the skyrmion bubble energy with electric field (Figure 3c,d) is calculated using the experimental 𝑀 𝑠 variations with electric field and a DW energy variation adjusted to fit the observed skyrmion bubble nucleation/annihilation rates (see supplementary S10). For positive voltage, skyrmion bubbles vanish due to a strong increase of the nucleation energy 𝐸 𝑛 combined with a decrease of the annihilation barrier:

nucleation is prevented and annihilation of bubbles is thermally activated. On the contrary, for a negative electric field, the nucleation energy 𝐸 𝑛 is strongly reduced whereas 𝐸 𝑎 becomes very large, thus favoring the creation of stable skyrmion bubbles. Consequently, the electric field can be used similarly to an applied magnetic field to switch skyrmion bubbles on and off by modifying the skyrmion stability. In our system the observed switch effect is associated to strong electric field variations of 𝑀 𝑠 and 𝜎 𝑊 only expected near 𝑇 𝐶 . However, the variation of 𝑀 𝑠 and 𝜎 𝑊 have opposite effects: an increase of 𝑀 𝑠 /𝜎 𝑊 is decreasing/increasing 𝐸 𝑛 . Consequently, the skyrmion switch effect should be much more efficient at fixed 𝑀 𝑠 . To check this we calculated in Figure 3f the skyrmion bubble energy in the case of a thicker 𝐶𝑜 layer of 0.6 nm, which is expected to show negligible 𝑀 𝑆 variation with electric field (higher 𝑇 𝐶 ) but a significant 𝐾 𝑒𝑓𝑓 change [START_REF] Bernand-Mantel | Electric-Field Control of Domain Wall Nucleation and Pinning in a Metallic Ferromagnet[END_REF] . We see that the skyrmion switch effect (i.e. reducing the nucleation barrier and increasing the annihilation barrier) can be efficient at fixed 𝑀 𝑠 with a 𝜎 𝑊 variation as small as 5%.

In this case, the skyrmion bubbles are scaled down to a few 100 nm of diameter. Below this size, our thin wall analytical model is no more valid and further simulations are necessary to study the skyrmion swich effect on a nanometer sized skyrmion in a confined geometry.

In conclusion, we have reported the presence at RT of micron sized skyrmion bubbles in a 𝑃𝑡/𝐶𝑜/𝑜𝑥𝑖𝑑𝑒 trilayer. The DMI strength, which is of the order of 1-2 mJ/m², reduces the DW energy by a factor of ~5 with respect to the classical Bloch case and enables the stabilization of skyrmion bubbles. The observed skyrmion bubbles density variation along the wedge and versus electric field is well described by an analytical isolated bubble model. We demonstrate the efficient and reproducible electric field writing and deleting of skyrmion bubbles. These results constitute a potentially important milestone towards the use of skyrmions for memory or logic devices.

Methods

Sample preparation

The 𝑃𝑡(3 nm)/𝐶𝑜 (0.6 nm)/𝐴𝑙 0.76 nm to 𝐴𝑙 1.76 nm film was deposited by magnetron sputtering on a Si wafer with a 500 nm thick thermal 𝑆𝑖𝑂 2 top layer. The 𝐴𝑙 layer was oxidized in-situ by a 𝑂 2 plasma.

A ~90 nm 𝐻𝑓𝑂 2 layer is deposited ex-situ by Atomic Layer Deposition. Indium Tin oxide electrodes deposited by sputtering was patterned by lift-off in the form of 50 × 800 μm rectangles.

Analytical model of the labyrinthine domains width

The theoretical labyrinthine domain width 𝐿 has been calculated using a parallel band domains model in the case where the DW width is negligible compared to the domain diameter. For ultra-thin films, when the thickness is much smaller than the characteristic dipolar length [START_REF] Kaplan | The Domain Structure in Ultrathin Magnetic Films[END_REF][START_REF] Kooy | Experimental and Theoretical Study of the Domain Configuration in Thin Layers of BaFeO[END_REF] : 𝑡 ≪ 𝐿 0 with 𝐿 0 = 𝜎 𝑊 /𝜇 0 . 𝑀 𝑠 2 where 𝑀 𝑆 is the saturation magnetization and 𝜎 𝑊 the domain wall energy, the domain width 𝐿 can be written as 𝐿 = 𝛼. 𝑡. exp(𝜋𝐿 0 /𝑡) where 𝑡 is the magnetic film thickness and 𝛼 a constant 𝛼 = exp ( 𝜋𝑏 2 + 1)~0.955 where 𝑏 is the numerical evaluation of a series [START_REF] Kaplan | The Domain Structure in Ultrathin Magnetic Films[END_REF] . Using this expression we obtain the domain wall energy 𝜎 𝑊 = ln (𝐿/(𝛼. 𝑡)). 𝜇 0 . 𝑀 𝑠 2 .t/ 𝜋.

Analytical "thin wall" bubble model

We have used an analytical formula to approximate the energy of an individual bubble of radius 𝑅 relative to the saturated state [START_REF] Bobeck | Properties and Device Applications of Magnetic Domains in Orthoferrites[END_REF][START_REF] Malozemoff | Magnetic Domain Walls in Bubble Materials[END_REF] . The model used here is valid for 𝑄 = 𝐾 𝑢 /𝐾 𝑑 > 1, where 𝐾 𝑢 is the uniaxial (surface) out of plane anisotropy and 𝐾 𝑑 the magnetic dipolar energy. In our case we have 𝑄 = 0.83/0.54 = 1.53. It was shown that this standard theory is well applicable for 𝑄 > 1.5. The energy of an individual bubble relative to the saturated state can be written as:

∆𝐸 𝑠𝑏 = 𝜎 𝑤 . 𝑡. 2. 𝜋. 𝑅 + 2𝜇 0 𝑀 𝑠 . 𝐻. 𝑡. 𝜋. 𝑅 2 -𝜇 0 𝑀 𝑠 2 . 𝜋. 𝑡 3 . 𝐼(𝑑)

where 𝜎 𝑤 is the DW energy, 𝑡 the layer thickness, 𝑅 the bubble radius, 𝑀 𝑠 the saturation magnetisation, 𝜇 0 𝐻 the applied magnetic field, and 𝐼(𝑑) defined as follow:

𝐼(𝑑) = - 2 3𝜋 𝑑[𝑑 2 + (1 -𝑑 2 )𝐸(𝑢 2 )/𝑢 -𝐾(𝑢 2 )/𝑢]
where 𝑑 = 2𝑅/𝑡, 𝑢 2 = 𝑑 2 /(1 + 𝑑²) and 𝐸(𝑢) and 𝐾(𝑢) are elliptic integrals defined as:

𝐸(𝑢) = ∫

  . The calculated nucleation and annihilation energies are shown in Figure1d,e. For this calculation, we have used the 𝑀 𝑠 , 𝑡 and 𝜎 𝑊 values and their variations along the wedge estimated experimentally using a combination of magneto-optical Kerr and VSM-SQUID magnetometry (see supplementary S1-S3). As our characterization setup was adapted to study the 100 ms to few s timescale, we draw as a guide to the eye, a white line indicating the RT thermal activation frontier, 𝐸/𝑘 𝐵 𝑇 293𝐾 = 23 (value at witch nucleation events are occurring at a rate  > 1 Hz). The analytical skyrmion bubble model reproduces several features observed in Figure1a. In Figure1d, for a constant magnetic field, the increase of the nucleation barrier 𝐸 𝑛 with the DW energy explains the decrease in the bubble density observed in Figure1aas the 𝐶𝑜 gets thicker. In the zone indicated by the number 2 in Figure1d, skyrmion bubbles are observed despite the increase of 𝐸 𝑛 above the thermal activation frontier. Those skyrmion bubbles are metastable skyrmions bubbles nucleated at lower magnetic field or stabilized by local magnetic inhomogeneities, which might locally lower the nucleation barrier. In the position indicated by the number 1 in Figure1d,e, both nucleation and annihilation are expected with a rate  > 1 Hz. This thermal instability of the skyrmion bubbles is observed in a time dependent measurement (see supplementary S7). On the contrary, the skyrmion bubbles present a high thermal stability (𝐸 𝑎 /𝑘𝑇 > 23) in a large part of the diagram, below the white line in Figure1e. In this region, no annihilation of skyrmion bubbles are observed from a few minutes and up to much longer timescales. The position 𝑝 1 near this white line has been selected for the observation of skyrmion bubbles under electric field. The behavior of skyrmion bubbles under electric field gating is studied through a transparent 𝐼𝑇𝑂 electrode (see Figure1b). Polar Kerr images recorded for different applied electric fields are presented in Fig 2 a-c. A distribution of bubble diameters due to the presence of local inhomogeneities and pinning

Figure 1 :

 1 Figure 1: (a) Polar Kerr magnetic images of the labyrinthine and skyrmion bubble states in the Pt/Co/oxide trilayer. The image is made of 10 images of 30 x 550 µm. Each image is recorded a few seconds after applying an out of plane magnetic field. The magnetic field was varied from 0 to 0.3 mT. The Co thickness increases from left to right. (b) Schematic representation of the device: the Pt/Co/oxide trilayer is covered by a 90 nm HfO2 layer and a top Indium Tin Oxide (ITO) electrode. The Co thicknesses variation is of the order of 1% over the 500 µm wide observed region. (c) Characteristic domain width 𝐿 as function of the wedge position for zero applied magnetic field extracted by fast Fourier transform on polar Kerr images. (d) and (e) Skyrmion bubble nucleation (d) and annihilation energies (e) calculated with the isolated bubble model as function of out of plane magnetic field and DW energy using the 𝐻, 𝑡, 𝑀 𝑠 and 𝜎 𝑤 parameters estimated experimentally.

Figure 2 :Figure 3 :

 23 Figure 2: (a) to (c) Polar Kerr magnetic images of the electric field control of skyrmion bubbles density in the Pt/Co/oxide trilayer under a static 0.15 mT perpendicular magnetic field, recorded through the transparent ITO electrode near position 𝑝 1 with different applied voltages. The sample is first saturated with a higher magnetic field, then the magnetic field is fixed to 0.15 mT and the electric field is varied from -5 V to +10 V. (d) and (e) Differential images obtained from (b)-(a)=(d) and (c)-(b)=(e). The blue/red bubbles correspond to objects which appeared/disappeared during the few seconds separating the images acquisitions. (f) Number of skyrmion bubbles extracted from images with twice the size of (a),(b) and (c). (g) Analytical calculation of the energy difference between a saturated state and a single isolated magnetic bubble as function of the bubble diameter with the parameters 𝑡 = 0.468 nm, 𝜇 0 𝐻 = 0.15 mT, 𝑀 𝑠 = 0.92+/-0.05 MA/m and 𝜎 𝑤 = 1.33+/-0.16 mJ/m² corresponding to respectively 0 V, +5 V and -5 V. The nucleation 𝐸 𝑛 and annihilation 𝐸 𝑎 energies are indiczted for the +5 V case. The 𝑡, 𝐻, 𝑀 𝑠 and 𝜎 𝑤 parameters used for the simulation are determined experimentally (see supplementary S3, S8)

Figure S1 :

 S1 Figure S1 : Schematic representation of the sample A1 with the 0-5 and 𝒑 𝟎 -𝒑 𝟑 positions used in the text.

Figure S2 :

 S2 Figure S2 : VSM-SQUID measurements of the remanent magnetization of a sample similar to A1 as function of the temperature. Each point is the remanent magnetization extracted from M(H) loops.

Figure S3 : a ,

 a Figure S3 : a, Polar Kerr loop maximum amplitude recorded every 160 µm along the wedge starting from a position around 𝒑 𝟎 divided by the Kerr loop maximum amplitude at position 𝒑 𝟎 . The corresponding magnetization values are indicated on the right y axis. b, Polar Kerr loops versus in-plane magnetic field recorded at positions 1 (black squares) and 3 (blue triangles) indicated on a. SW fit (black and blue lines, see S4).

Figure S4 :

 S4 Figure S4 : Experimental values of the DW energy along the wedge calculated using the experimental 𝑴 𝐬 , 𝒕 and 𝑳 values (black squares) and polynomial fit (red line) (see S3).

Figure S5 :

 S5 Figure S5: Definition of the angles between the magnetization, the applied magnetic field, the film plane and the observation direction (represented by an eye).

Figure S6 :

 S6 Figure S6: (a) Polar Kerr image of the labyrinthine domains under electric field at zero applied magnetic field.The dashed rectangle correspond to the zone where the Kerr signal of (b) was recorded b, Normalized polar Kerr signal for in plane applied magnetic field at positions 𝑝 1 on the wedge under applied voltages of 0 V (black circles) and -10 V (blue squares) superimposed with the SW fit (blue and black lines). c, Experimental (squares and dots) DW energies as function of the wedge position and for two applied voltages calculated from the analytical expression of the labyrinthine domain width using the 𝐿 values extracted by Fast Fourier transform from a. The black and blue lines are linear fits.

Figure S7 :

 S7 Figure S7: Normalized polar Kerr signal for in plane applied magnetic field at positions 𝑝 1 under applied voltages of +20 V (black circles) and -20 V (blue circles) superimposed with the SW fits (black and blue lines).

Figure S8 :

 S8 Figure S8: Isolated skyrmion bubble energy as function of the bubble diameter calculated with 𝑡 = 0.468 nm, μ 0 𝐻 = 0.15 mT and the magnetization values corresponding to (a) the +20 V case with 𝑀 𝑠 = 1.07 MA/m and (b) the -20 V case with 𝑀 𝑠 = 0.57 MA/m. The DW energy parameters are indicated in the graphs.

  2 , the saturation magnetization 𝑀 𝑆 , the DW energy 𝜎 𝑊 and 𝛼 a numerical constant (see Methods). The experimental value of the DW energy 𝜎 𝑊 can be calculated from this equation with the measured values of 𝐿, t and 𝑀 𝑆 (see supplementary S1-S3). At position 𝑝 0 , 𝐿 = 2.9 μm, 𝑡 = 0.47 nm, 𝑀 𝑆 = 0.93 MA/m and we find 𝜎 𝑊 = 1.44 mJ/m². This value can be compared to the classical Bloch wall energy 𝜎 𝑤 𝐵𝑙𝑜𝑐ℎ = 4 √ 𝐴 𝑒𝑥 𝐾 𝑒𝑓𝑓 .

  𝑡, 𝑀 𝑠 , 𝜎 𝑤 ) = 𝜎 𝑤 . 𝑡. 2. 𝜋. 𝑅 + 2𝜇 0 𝑀 𝑠 . 𝐻. 𝑡. 𝜋. 𝑅 2 -𝜇 0 𝑀 𝑠 2 . 𝜋. 𝑡3 . 𝐼(𝑑)

Acknowledgement

This work was supported by the French National Research Agency (ANR) under the project ELECSPIN ANR-16-CE24-0018. We acknowledge the support from the Nanofab facility and the pôles magnetometrie and ingénierie experimentale from Institut Néel. We thank L. Cagnon for ALD deposition. The authors thanks A. Thiaville, O. Fruchard, M. Chshiev, O. Boulle, I. M. Miron, N. Reyren, V. Cros, A. N. Bogdanov and B. A. Ivanov for fruitful discussions.

Supporting information

Sample characterisation. Current induced skyrmion bubbles motion. Real time room temperature spontaneous nucleation and annihilation of skyrmion bubbles. Sample characterisation under electric field. Switching of skyrmion bubbles with electric field. Estimation of the electric field variation of the domain wall energy using the isolated bubble model.

Competing financial interests

The authors declare no competing financial interests.

The skyrmion switch: turning magnetic skyrmion bubbles on and off with an electric field Marine Schott 1,2,3,4 , Anne Bernand-Mantel 1,2 *, Laurent Ranno 1,2 , Stéfania Pizzini 1,2 , Jan Vogel 1,2 , Hélène Béa 1,3,4 , Claire Baraduc 1,3,4 , Stéphane Auffret 1,3,4 , Gilles Gaudin 1,3,4 and Dominique Givord 1,2

1 Univ. Grenoble Alpes, F-38000 Grenoble, France 2 CNRS, Inst. NEEL, F-38000 Grenoble, France 3 CNRS, SPINTEC, F-38000 Grenoble, France 4 CEA, INAC-SPINTEC, F-38000 Grenoble, France *anne.bernand-mantel@neel.cnrs.fr

S1. VSM-SQUID magnetization measurements.

The 𝑃𝑡(3 nm)/𝐶𝑜 (0.6 nm)/𝐴𝑙 0.76 nm to 𝐴𝑙 1.76 nm + 𝑂 2 plasma sample has been cut to carry out VSM-SQUID measurements. The piece where the labyrinthine domains discussed in the paper appear is called sample A1 and has dimensions of ~1 mm × 4 mm. The A1 sample is represented schematically in Fig. S1 where the 𝑝 0 , 𝑝 1 , and 𝑝 3 positons are the one indicated also in Fig. 1a. The magnetic moment

is measured with a Quantum Design MPMS 3 VSM-SQUID equipment. The sample A1 was measured 4 times at room temperature and we find 𝑀 𝑠 . 𝑡 = 4.2 ± 0.4 × 10 -4 A/m.m. As a part of the 𝐶𝑜 has been oxidized by 𝑂 2 plasma in the sample A1, the 𝐶𝑜 thickness is lower than the initial thickness 𝑡 = 0.6 nm.

The temperature dependent magnetization measurements, carried out in the 100-360K temperature range, on a sample similar to A1, are presented in Fig. S2. Each point correspond to the remanent magnetisation extracted from a M(H) loop. A Curie temperature 𝑇 𝐶 = 325 K is found using a combination of two power laws :

S2. Polar Kerr measurements along sample A1.

To estimate the 𝑀 𝑠 variation along the sample A1, local Kerr measurements were carried out every 160 µm. In Fig. S3b, we show the polar Kerr measurements with in plane applied magnetic field. This measurement correspond to the out of plane projection of the magnetization while it rotates reversibly outside its easy axis toward its hard axis direction (see Fig. S5a). We define the maximum Kerr signal amplitude as the difference between the maximum and the minimum intensities. This Kerr signal amplitude is proportional to the Kerr rotation, itself proportional to 𝑀 𝑠 . 𝑡 𝐶𝑜 . The Kerr signal maximum amplitude is plotted in Fig. S3a The thickness is fixed also to its lower bound 𝑡 = 0.47 nm. To estimate the variation of the Co thickness along the wedge we made the assumption that the Co thickness variation due to oxidation is a linear projection of the aluminium thickness variation (see Fig S1). The deposition is done by sputtering in a system with square targets. The wedge is obtained by progressively moving the shutter during the deposition. To calibrate the thickness variation over the 10 cm width of the substrate, a thicker wedged layer has been deposited and the thickness has been measured regularly using X-ray reflectivity along the sample. Assuming that the slope varies proportionally with the thickness, we deduce a slope of 0.1 nm/cm for the aluminium wedge, and consequently the Co layer in the overoxidized zone, assuming that the oxygen plasma is homogeneous.

In Fig. S4, we show the domain wall energy and its variations, calculated from the analytical expression of the labyrinthine domains width defined in the Methods : 𝜎 𝑤 = ln(𝐿/(0.955. 𝑡)). 𝜇 0 𝑀 𝑠 ². t/π using the 𝑀 𝑠 and 𝑡 values given previously and the experimental domain width 𝐿 values extracted by fast Fourier transform from Fig. 1a which is plotted in Fig. 1c. These values are fitted with a polynomial expression:

𝜎 𝑤 = 1.44. 10 -3 -4.53. 10 -7 𝑥-5.43. 10 -10 𝑥² where 𝑥 is the position in µm (Fig. S4). The table T1 summarizes the experimental 𝑀 𝑠 , 𝑡 and 𝜎 𝑤 values defined here and their variations along the wedge.

These values will be used as an input for the isolated skyrmion bubble simulations in Fig. 1 d ande. We carry out simulations of the 𝑀(𝐻) loops in the hard axis direction based on the Stoner-Wohlfarth (SW) model which allows to reproduce the magnetization rotation and estimate 𝐾 𝑒𝑓𝑓 /𝑀 𝑠 . We assume that the system has a uniaxial easy axis perpendicular to the film surface (Fig. S5). The energy per unit volume is: : 𝐸 = 𝐾 𝑒𝑓𝑓 𝑠𝑖𝑛 2 (𝜃) -𝜇 0 𝐻𝑀 𝑠 𝑐𝑜𝑠(𝜃 -𝜑) where 𝐾 𝑒𝑓𝑓 is the uniaxial anisotropy in 𝐽/𝑚 3 , 𝜃 the angle between the magnetization and the easy axis and 𝜑 the angle between the magnetic field and the easy axis (see Fig. S5). An angle 𝛼 is introduced to take into account the misalignment between the applied magnetic field, the film plane and the direction perpendicular to the observation direction (Fig. S5b). In Fig. S3b

S5. Real time room temperature spontaneous demagnetization at zero applied magnetic field.

The video V1 displays in real time the polar Kerr images of the room temperature spontaneous demagnetization. The sample was saturated with an out of plane magnetic field which is turned to zero just after the beginning of the video. We observe a spontaneous demagnetization of the sample in few seconds via the nucleation and growth of magnetic domains. When the sample is demagnetized, the thermally activated movement of the DW is still visible. Such thermal melting of stripe domains has been observed in systems near the spin reorientation transition [3][4][5][6][7][8][9] .

S6. Current induced skyrmion bubbles motion.

The sample was etched to form a Hall cross in order to check the unidirectional skyrmion bubbles displacement with electrical current. The channels of the Hall cross are 10 µm large. A DC current of 0.15 mA is applied between the top right and bottom contact. The video V2 displays in real time the polar Kerr image of the Hall cross at room temperature under a ~0.15 mT constant magnetic field applied in the black direction. The skyrmion bubbles have a steady motion against the electron flow for a current density of 5 × 10 9 𝐴/𝑚². Some bubbles get pinned and/or distort into elongated domains reminiscent of stripe domains. This result is similar to what has been obtained in other systems presenting skyrmion bubbles [10][11][12] . A more detailed bubble motion study was not carried out in our work.

S7. Real time room temperature spontaneous nucleation and annihilation of skyrmion bubbles.

In video V3, the room temperature real time polar Kerr image near the number 1 indicated in Fig. 1d-e under a constant out of plane magnetic field of ~0.25 mT applied in opposite direction to the bubbles magnetization. Spontaneous nucleation and annihilation of white skyrmion bubbles with a rate > 1 Hz are observed.

S8. Experimental determination of the electric field variation of 𝑴 𝒔 , 𝑲 𝒆𝒇𝒇 and 𝝈 𝑾 using Kerr measurements and the Stoner-Wohlfarth model.

Polar Kerr loops with in plane applied magnetic field have been recorded at position 𝑝 1 = 215 μm under electric field to characterize the variations of 𝑀 𝑠 and 𝐾 𝑒𝑓𝑓 . The 𝑀 𝑠 value at position 𝑝 1 , determined in paragraph S3 is 𝑀 𝑠 = 0.92 MA/m. We use the SW model described in the paragraph S4 to deduce 𝐾 𝑒𝑓𝑓 and extract the 𝑀 𝑠 and 𝐾 𝑒𝑓𝑓 variations with electric field. In Fig. S6b, the normalised polar Kerr signal for in plane applied magnetic field measured at positions 𝑝 1 with 0 V and -10 V applied voltages are superimposed with the SW fits obtained with the parameters 𝑀 𝑠 = 0.92 MA/m (0 V), 𝑀 𝑠 = 0.82 MA/m (-10 V), 𝐾 𝑒𝑓𝑓 = 0.27 MJ/m 3 (0 V), 𝐾 𝑒𝑓𝑓 = 0.22 MJ/m 3 (-10 V) , 𝜑 = 84° and , 𝛼 = 5°. In Fig. S6a, a polar Kerr image is showing the effect of the electric field on the labyrinthine domains: the domains size is shifted by ~ 400 μm toward the thinner wedge direction for negative voltage (electron depletion in Co). This observation is consistent with the equivalence between a reduction of the electron density and an oxidation increase observed in our previous work 13 . The labyrinthine domain width 𝐿 and its electric field variation is extracted from the image of Fig. S6a

S9. Switching skyrmion bubbles with electric field.

The video V4 displays, in real time, at room temperature, the polar Kerr image of the skyrmion switch effect at position 𝑝 1 . The sample was saturated under a high out of plane magnetic field in the opposite direction (black) to the bubbles magnetization before the beginning of the video, then a 0.15 mT out of plane magnetic field is applied in the black direction and kept constant. During the video the electric field is switched every two seconds between -20 V and + 20 V starting with + 20 V.

S10. Estimation of the electric field variation of the DW energy 𝝈 𝑾 for high electric field using the isolated bubble model.

At high electric fields the labyrinthine domains disappear as the thermal activation is either too high (at -20 V), leading to a complete melting of the labyrinthine domains, or too low (at +20 V), preventing the system to demagnetize and the labyrinthine domains to form. Consequently, the experimental determination of the DW energy 𝜎 𝑤 using the labyrinthine domain width 𝐿 was not possible for high electric fields. To estimate 𝜎 𝑤 for these high voltages, we applied the isolated bubble model, using as an input the 𝑀 𝑠 values for high electric field obtained in paragraph S8 from Kerr measurements under electric field. We start by the +20 V case. At this voltage all the skyrmion bubbles disappear, as observed in Fig. 3b and video V4. This means that the annihilation energy barrier (defined in Fig. 2g) is smaller than the RT thermal energy limit 𝐸 𝑎 /𝑘 𝐵 𝑇 293𝐾 < 23, and much smaller than the nucleation energy 𝐸 𝑎 ≪ 𝐸 𝑛 . This is true for DW energies 𝜎 𝑤 > 1.85 mJ/m² as shown in Fig. S8a. This gives a lower bound for the increase of 𝜎 𝑤 under electric field ∆𝜎 𝑤 𝑝𝑜𝑠 =𝜎 𝑤 (20 V)-𝜎 𝑤 (0 V)=1.85-1.33>0.52 mJ/m². For -20 V, we observe a fast nucleation of skyrmion bubbles in Fig. 3a and video V4 and we deduce a nucleation energy 𝐸 𝑛 /𝑘 𝐵 𝑇 293𝐾 < 23, an annihilation energy 𝐸 𝑎 /𝑘 𝐵 𝑇 293𝐾 > 23. This gives a DW energy 𝜎 𝑤 < 0.47 mJ/m² as shown in Fig. S8b and a lower bound for the decrease of 𝜎 𝑤 under electric field ∆𝜎 𝑤 𝑛𝑒𝑔 = 𝜎 𝑤 (0 V)-𝜎 𝑤 (-20 V)=1.33-0.47>0.86 mJ/m².