Michael Urban

Coupling of hydrodynamics and quasiparticle motion in collective modes of superfluid trapped Fermi gases

Keywords: numbers: 03.75.Ss, 03.75.Kk, 67.40.Bz

Coupling of hydrodynamics and quasiparticle motion in collective modes of superfluid trapped Fermi gases

M.

I. INTRODUCTION

Most of the current experiments involving trapped atomic Fermi gases focus on the BEC-BCS crossover. By changing the magnetic field around a Feshbach resonance, the scattering length a of the atoms can be varied from small positive values through very large values near the resonance to small negative values. For a > 0, k F a ≪ 1 (where k F denotes the Fermi momentum) the system can be considered as a Bose-Einstein condensate (BEC) of diatomic molecules. The crossover region, k F |a| > ∼ 1, is not yet very well understood from a theoretical point of view. Finally, on the other side of the resonance, when a < 0, k F |a| ≪ 1, the system should be in the BCS phase if the temperature is sufficiently low. However, the BCS critical temperature T c is extremely low, and very soon the magnetic field reaches the point where T c becomes smaller than the actual temperature T , and the system undergoes the phase transition to the normal (non-superfluid) phase.

One possibility to study the crossover experimentally is to measure the properties of certain collective oscillations. For example, the radial and axial breathing modes of a cigar-shaped trapped Fermi gas have been measured over the whole crossover region [1,2]. In these experiments one can observe how the frequencies and damping rates of the modes change from what one expects for a BEC to what one expects for a collisionless normal Fermi gas. Assuming that, except in the collisionless normal phase, hydrodynamics is valid, the measured frequencies can give some information on the equation of state in the crossover region.

However, this schematic picture is not completely accurate. Since the system is in a trap, there is no sharp transition from the superfluid to the normal phase. This can be seen as follows: The BCS critical temperature T c depends on the atom density ρ, and the density depends on the position r. In the center of the trap, the density ρ(r) and hence the local critical temperature T c (r) are higher than in the outer part of the trap. As a consequence, for a given temperature, the outer part gets already normal at a magnetic field where the inner part is still superfluid. To be more precise, a system in the BCS phase at finite temperature behaves effectively like a mixture of superfluid and normal components with densities ρ s and ρ n , respectively, which become ρ s = ρ, ρ n = 0 in the limit T = 0 and ρ s = 0, ρ n = ρ in the limit T ≥ T c . As a consequence, if 0 < T < T c (r = 0), the superfluid inner part of the trap behaves like a mixture of normal and superfluid components, while only the outer part with T c (r) < T is completely normal [3].

If the collision rate was high enough, also the normal component of the gas would behave hydrodynamically. Such a system could be described by Landau's twofluid hydrodynamics which has been applied to collective modes in trapped superfluid gases at finite temperature [4]. However, although in the recent experiments the transition to the normal phase seemed to occur at a value of k F |a| ≈ 2 [1] (i.e., the BCS phase has not really been reached), the system behaved already like a collisionless normal Fermi gas. Hence it seems to be clear that the normal component cannot be treated in terms of hydrodynamics, but a description in terms of a Vlasov equation is required.

We note that there are other approaches to the description of the collective modes at finite temperature. In particular, let us mention the quasiparticle random phase approximation (QRPA) [5,6], which can be seen as the linearized form of the time-dependent Bogoliubov-de Gennes (BdG) equations. However, for practical reasons this method is limited to systems with spherical symmetry and numbers of particles up to a few times 10 4 . Another disadvantage of this method is that it does not allow to include a collision term.

For the case of clean superconductors, a semiclassical transport theory taking the coupling between normal and superconducting components into account has been developed by Betbeder-Matibet and Nozières [7]. Transport theories of this type have also been used for describing the dynamics of superfluid 3 He [8,[START_REF] Vollhardt | The Superfluid Phases of Helium 3[END_REF]. In a preceding paper [10], we derived the semiclassical transport equations for the case of trapped atomic Fermi gases and applied them to the quadrupole mode of a gas in a spherical trap. We found that the presence of the normal component leads to a strong damping of the hydrodynamic collective mode. The same mechanism might explain the strong damping observed experimentally near the transition to the collisionless behavior [1]. However, in Ref. [10] we had to replace the gap ∆(r) by a constant in order to find an analytical solution of the transport equations. Due to this simplification, which cannot really be justified, the damping of the hydrodynamic mode at a given temperature was much weaker than that obtained in QRPA calculations [6].

In the present paper we will work out a numerical method which allows us to treat the realistic rdependence of the gap. In addition, the method is very versatile and allows to treat much more general cases than can be solved analytically in the constant-gap approximation. The basic idea is to replace the continuous phase-space distribution function of the quasiparticles by a sum of a finite number of delta functions in phase space, called "test particles." In the normal phase, the testparticle method is routinely used for solving the Vlasov equation, e.g. for simulating heavy-ion collisions in nuclear physics [11]. It has also been applied to the simulation of the dynamics of normal trapped atomic Fermi gases with collision term [12] and of Bose-Fermi mixtures [13]. However, to our knowledge, the test-particle method has not yet been used in the context of superfluid systems, and in fact the numerical difficulties are quite different from those encountered in the usual applications.

The article is organized as follows. In Sec. II, we give a brief summary of the transport equations for the BCS phase and their linearization in the case of small deviations from equilibrium. We also give arguments why some terms which appear in the equations can be neglected. In Sec. III we introduce the test-particle method for the case of small oscillations around equilibrium. We describe in detail a number of tricky points we encountered during the implementation of the method, in particular the calculation of the test-particle trajectories, the generation of the test-particle distribution in phase space, and the initialization after a delta-like perturbation. In Sec. IV we present the first results obtained with the help of this method, again for the quadrupole mode in a spherical system. Finally, in Sec. V, we summarize and draw our conclusions.

II. TRANSPORT EQUATIONS FOR THE BCS PHASE A. Summary of the kinetic equations

In this subsection we will give a brief summary of the kinetic equation approach developed by Betbeder-Matibet and Nozières [7] for the case of clean superconductors and adapted to the case of trapped atomic Fermi gases in Ref. [10]. In this paper we will only give the final equations. For details of the derivations, see Ref. [10] and the footnote [27].

We consider a dilute gas of fermionic atoms of mass m in two equally populated hyperfine states ↑ and ↓, trapped by an external potential V ext and interacting via an attractive short-range interaction which leads to a scattering length a < 0. The corresponding classical mean-field hamiltonian (minus the chemical potential µ) reads

h(r, p) = p 2 2m + V (r) -µ , (1)
where V denotes the sum of the external and the Hartree potential,

V (r) = V ext (r) + V Hartree (r) = V ext (r) + gρ(r) . (2)
In the latter equation, g = 4π 2 a/m denotes the coupling constant and ρ is the density per spin state. The Vlasov equation (without collision term) for the phasespace distribution function ̺(r, p) in the normal phase can be written in the compact form

̺ = {h, ̺} , (3)
where {•, •} denotes the Poisson bracket. One way to derive this equation is to perform a Wigner-Kirkwood expansion up to order of the time-dependent Hartree-Fock equation [11,[START_REF] Ring | The Nuclear Many-Body Problem[END_REF].

In the superfluid phase the derivation of an analogous transport equation is much more complicated due to the presence of the complex order parameter (gap) ∆(r) whose phase describes the collective motion of the Cooper pairs. In addition to the density matrix ̺, there exists now an anomalous density matrix (pairing tensor) κ. The gap ∆ and the anomalous density are related by the gap equation

∆(r) = -g d 3 p (2π ) 3 κ(r, p) - ∆(r) p 2 /m . ( 4
)
The time-dependence of ̺ and κ is governed by the timedependent Hartree-Fock-Bogoliubov or BdG equations. As in the normal phase, the semiclassical transport theory can be derived from these equations by performing a Wigner-Kirkwood expansion up to order . However, it turns out that it is necessary to introduce a gauge transformation with a phase φ(r) that makes the order parameter ∆ real. This corresponds to a transformation into the local rest frame of the Cooper pairs, which are moving with the collective velocity v coll (r) = -( /m)∇φ(r). The effect of this transformation is to change the gap ∆, the single-particle hamiltonian h, the normal and anomalous density matrices ̺ and κ according to

∆(r) = ∆(r)e 2iφ(r) ≡ |∆(r)| , (5) h(r, p) = h[r, p -∇φ(r)] -φ(r) , (6)
̺(r, p) = ̺[r, p -∇φ(r)] , (7)
κ(r, p) = κ(r, p)e 2iφ(r) .

Roughly speaking, the phase φ determines the dynamics of the superfluid component of the system, while the dynamics of the normal component, consisting of thermally excited quasiparticles, has to be described separately. The distribution of these quasiparticles, denoted by ν(r, p), obeys the following equation of motion:

ν = {E, ν} . (9)
This equation looks formally very similar to the Vlasov equation ( 3), except that the hamiltonian h is replaced by the quasiparticle energie E, which is defined as

E = h2 ev + ∆2 + hod . ( 10
)
Throughout this article, the indices "ev " and "od" denote the time-even and time-odd parts of a phase-space function, i.e., the parts which are even and odd in p, respectively. The quasiparticle-distribution function ν is related to the normal and anomalous density matrices in the new gauge, ̺ and κ, by

̺ = 1 2 - hev 2E ev (1 -2ν ev ) + ν od , (11)
Re κ = ∆ 2E ev (1 -2ν ev ) . (12)
The Vlasov-like equation [START_REF] Vollhardt | The Superfluid Phases of Helium 3[END_REF] has to be complemented with an equation of motion for the phase φ. It turns out that φ has to be determined from the continuity equation

ρ(r) + ∇ • j(r) = 0 , (13)
where the density ρ and the current j are given by

ρ(r) = d 3 p (2π ) 3 ̺(r, p) , (14)
j(r) = d 3 p (2π ) 3 p m ̺(r, p) -m ρ(r)∇φ(r) . (15)

B. Linearization around equilibrium

Let us now assume that the external potential V ext can be written as

V ext = V 0ext + V 1ext , (16)
where V 0ext is time-independent and V 1ext is a small perturbation. The equilibrium quantities (corresponding to the potential V 0ext ) will be marked by an index "0". In particular, since in equilibrium the gap can be chosen to be real, we have

φ 0 = 0 , h0 = h 0 , ∆0 = ∆ 0 . (17)
The quasiparticle distribution function in equilibrium is given by

ν 0 (r, p) = f [E 0 (r, p)] , (18)
where f (E) denotes the Fermi function,

f (E) = 1 e E/(kB T ) + 1 . ( 19
)
Our aim is to calculate the small deviations from equilibrium induced by the perturbation V 1ext , which will be marked by an index "1". To that end we linearize the transport equation ( 9) for the quasiparticle-distribution function,

ν1 -{E 0 , ν 1 } = f ′ (E 0 ){E 1 , E 0 } , (20)
where f ′ (E 0 ) = df /dE 0 . The deviation of the quasiparticle energy, E 1 , which appears on the r.h.s., depends itself on ν 1 through the deviation of the Hartree field, gρ 1 , and the deviation of the gap, ∆1 . Expressing ρ 1 and ∆1 in terms of ν 1 , we can write Eq. ( 20) as

ν1 -{E 0 , ν 1 } = - f ′ (E 0 ) m -p • ∇ V 1ext + gρ 1ν -φ1 1 + gA + ∆ 0 E 2 0 p • ∇ ∆ 0 (V 1ext + gρ 1ν -φ1 ) 1 + gA + h 0 E 2 0 p • ∇ ∆ 0 ∆ 1ν gA + m h 0 E 0 (p • ∇) 2 φ 1 - h 0 E 0 (∇V 0 ) • ∇φ 1 - ∆ 0 E 0 (∇∆ 0 ) • ∇φ 1 . ( 21
)
where ρ 1ν and ∆ 1ν are the quasiparticle contributions to ρ 1 and ∆1 ,

ρ 1ν (r) = d 3 p (2π ) 3 h 0 (r, p) E 0 (r, p) ν 1 (r, p) , (22)
∆ 1ν (r) = g d 3 p (2π ) 3 ∆ 0 (r) E 0 (r, p) ν 1 (r, p) , (23)
while A(r) is a function which depends only on equilibrium quantities. The explicit expression for the function A(r) reads

A(r) = mp F (r) 2π 2 3 [1 -ϕ(r)] , (24)
where the local Fermi momentum p F (r) is defined as usual by p 2 F (r)/(2m) = ǫ F (r) = µ -V 0 (r), and the temperature dependence of A(r) is governed by the function

ϕ(r) = -dξ ξ 2 E 2 ξ f ′ (E ξ ) , (25)
with E ξ = ξ 2 + ∆ 2 0 (r). In the two limiting cases T = 0 and T ≥ T c (r), the function ϕ(r) takes the values 0 and 1, respectively. As a consequence, A(r) = 0 if T ≥ T c (r).

In order to determine the phase φ 1 , we also linearize the continuity equation (13):

ρ1 (r) + ∇ • j 1 (r) = 0 . (26)
Again, we express all quantities in terms of equilibrium quantities and the unknown quantities ν 1 and φ 1 . According to Eq. ( 15), the current j 1 can be decomposed into quasiparticle and superfluid contributions,

j 1 (r) = j 1ν (r) -m ρ 0 (r)∇φ 1 (r) = 0 , (27)
where the quasiparticle contribution is given by

j 1ν (r) = d 3 p (2π ) 3 p m ν 1 (r, p) (28)
(note that only the time-odd part of ν 1 contributes to the integral, and ν 1od = ̺1od ). The time derivative ν1 which appears when one writes down the explicit expression for ρ1 can be eliminated with the help of Eq. ( 21). As a result, the l.h.s. of the continuity equation becomes

ρ1 + ∇ • j 1 = A 1 + gA φ1 -V1ext - 2π 2 3 mp F + g m ∇ • (ρ 0 ∇φ 1 ) + g∇ • j 1ν + ∆ 0 A ∇ • d 3 p (2π ) 3 ∆ 0 E 2 0 p m ν 1 . ( 29
)
As noted in Ref. [10], the continuity equation is trivially satisfied in the normal phase (T ≥ T c ). This becomes evident if its l.h.s. is written in the form (29), since in the normal phase we have ∆ 0 = 0 and A = 0.

C. Identification of important and unimportant terms

Eq. ( 21) is still very complicated. In order to simplify the problem, let us look more closely at the different terms in order to see if some of them are less important than others. The basic assumption being that ∆, k B T c , and k B T are much smaller than ǫ F , the distribution function is sharply peaked near the Fermi surface. Under this condition it is useful to express the distribution function ν in terms of the variables r, ξ, and p instead of r and p, where

ξ = h 0 (p, r) ≈ v F (r)[|p| -p F (r)] , p = p |p| , (30)
and v F (r) = p F (r)/m. In terms of these variables, ν is sharply peaked near ξ = 0, and the relevant values of ξ are of the same order of magnitude as ∆, k B T c , and k B T . If ν 1 is written as a function of the new variables, the Poisson bracket on the l.h.s. of Eq. ( 21) becomes

{E 0 , ν 1 } = ∆ 0 E 0 v F p • (∇∆ 0 ) ∂ν 1 ∂ξ + ξ E 0 1 p F (∇V 0 ) • ∂ν 1 ∂p - ξ E 0 v F p • ∇ν 1 + ∆ 0 E 0 1 p F (∇∆ 0 ) • ∂ν 1 ∂p , (31)
with the short-hand notation

∂ν 1 ∂p = ∂ν 1 ∂ϑ p ∂p ∂ϑ p + 1 sin 2 ϑ p ∂ν 1 ∂ϕ p ∂p ∂ϕ p (32)
where ϑ p and ϕ p denote the angles characterizing the unit vector p = (sin ϑ p cos ϕ p , sin ϑ p sin ϕ p , cos ϑ p ).

In addition to the assumption ∆ ≪ ǫ F , our semiclassical theory requires that all quantities vary slowly in space, i.e., on a length scale L which should be larger than the coherence length v F /(π∆) [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF]. Then, using ∆ 0 ∼ E 0 ∼ ξ ∼ ∆, ∇ ∼ 1/L, ∂/∂ξ ∼ 1/∆, and ∂/∂p ∼ 1, all terms in Eq. ( 31) can be estimated to be of the order of magnitude (v F /L)ν 1 , except the last one, which is of the order (v F /L)(∆/ǫ F )ν 1 . Hence, the last term of Eq. ( 31) is negligible.

Let us now distinguish different kinds of contributions to ν 1 , depending on whether they are even or odd functions in ξ and p: ν 1oe : the part of ν 1 which is odd in ξ and even in p describes, roughly speaking, a change of the Fermi momentum, i.e., fluctuations of the density, and contributes to ρ 1ν ,

ρ 1ν ≈ mp F 2π 2 3 dΩ p 4π dξ ξ E ν 1oe , (33)
with dΩ p = sin ϑ p dϑ p dϕ p , while its contribution to ∆ 1ν is suppressed by one power of ∆/ǫ F and can be neglected.

ν 1eo : the part of ν 1 which is even in ξ and odd in p describes a shift of the Fermi sphere and therefore contributes to the current j 1ν ,

j 1ν ≈ p 2 F 2π 2 3 dΩ p 4π p dξ ν 1eo , (34)
and also to the other integral in the continuity equation (29).

ν 1ee : the part of ν 1 which is even in ξ and in p describes, roughly speaking, a local temperature fluctuation and leads to a non-vanishing value of ∆1 (via ∆ 1ν ),

∆1 ≈ - 1 1 -ϕ dΩ p 4π dξ ∆ 0 E ν 1ee , (35)
while its contribution to ρ 1ν is suppressed by one power of ∆/ǫ F and can be neglected.

ν 1oo : the part of ν 1 which is odd in ξ and odd in p gives only a negligible contribution to the current j 1ν (suppressed by one power of ∆/ǫ F ).

If one neglects the last term in Eq. ( 31), the Poisson bracket in Eq. ( 21) leads only to a coupling between ν 1eo and ν 1oe and between ν 1oo and ν 1ee . To be more specific, ν 1ee and ν 1oo do not contribute to the dynamics of ν 1oe and ν 1eo . Since we are interested in density oscillations and currents, which are determined by ν 1oe and ν 1eo , we might wonder if we could disregard completely ν 1ee and ν 1oo . To that end we have to check that also on the r.h.s. of Eq. ( 21) there is no term which couples the undesired quantities ν 1ee and ν 1oo to ν 1oe or ν 1eo . Actually, on the r.h.s. of Eq. ( 21) there is no term containing ν 1oo and only one term containing ν 1ee , namely the third one,

- f ′ (E 0 ) m h 0 E 2 0 p • ∇ ∆ 0 ∆ 1ν gA ≈ v F f ′ (E 0 ) ξ E 2 0 p • ∇(∆ 0 ∆1 ) .
(36) This term clearly contributes to ν1oo , but at least to leading order in ∆/ǫ F it does not contribute to ν1eo or ν1oe . In the continuity equation ( 29), ν 1ee and ν 1oo do not appear, i.e., the undesired quantities ν 1ee and ν 1oo do not contribute to the dynamics of φ 1 , either. We are therefore allowed to disregard them. Now, since we are not interested any more in ν 1ee and ν 1oo , we can remove all the terms on the r.h.s. of Eq. ( 21) which contribute only to the dynamics of these uninteresting quantities. As mentioned above, this is the case for the third term, Eq. ( 36), which contributes only to ν1oo . The last term on the r.h.s. of Eq. ( 21),

f ′ (E 0 ) m ∆ 0 E 0 (∇∆ 0 ) • ∇φ 1 (37)
can be omitted, too, since it contributes only to ν1ee .

In conclusion, we are left with a simplified version of Eq. ( 21), which reads

ν1 -{E 0 , ν 1 } = - f ′ (E 0 ) m -p • ∇ V 1ext + gρ 1ν -φ1 1 + gA + ∆ 0 E 2 0 p • ∇ ∆ 0 (V 1ext + gρ 1ν -φ1 ) 1 + gA + m h 0 E 0 (p • ∇) 2 φ 1 - h 0 E 0 (∇V 0 ) • ∇φ 1 . ( 38
)

III. TEST-PARTICLE METHOD A. Description of the method

The aim of the present work is to solve the Vlasov-like equation ( 9) for the quasiparticle-distribution function ν together with the continuity equation ( 13) for the phase of the order parameter with the help of the test-particle method, in analogy to the test-particle method which is used to solve the usual Boltzmann equation. The basic idea of this method is to replace the continuous distribution function ν(r, p) by a sum of delta functions in phase space,

ν(r, p; t) ∝ i δ[r -R i (t)]δ[p -P i (t)] , (39)
corresponding to a finite number of test particles, each of which follows the classical equation of motion

Ṙi = ∂E(R i , P i ; t) ∂P i , Ṗi = - ∂E(R i , P i ; t) ∂R i , (40)
as can be seen by inserting Eq. (39) into Eq. ( 9). Note that, contrary to the usual test-particle method, our test particles here cannot be identified with real particles but rather with Bogoliubov quasiparticles. In its general form, the test-particle method can be applied to situations far from equilibrium. However, here we are only interested in the linear-response regime, i.e., in the limit of small deviations from equilibrium. In this case it is possible to formulate the method in such a way that only the classical trajectories corresponding to the unperturbed system appear.

To that end, we make the following ansatz for the deviation of the distribution function from equilibrium:

ν 1 (r, p; t) = -y(r, p; t)f ′ [E 0 (r, p)] . (41)
Inserting this into the linearized transport equation (38), we obtain the following equation of motion for the function y:

ẏ(r, p; t) -{E 0 (r, p), y(r, p; t)} = F (r, p; t) , (42)
where

F (r, p; t) = - p m • ∇ V 1ext + gρ 1ν -φ1 1 + gA + ∆ 0 E 2 0 p m • ∇ ∆ 0 (V 1ext + gρ 1ν -φ1 ) 1 + gA + h 0 E 0 p m • ∇ 2 φ 1 -m h 0 E 0 (∇V 0 ) • ∇φ 1 . (43)
Denoting by R(r, p; t) and P(r, p; t) the classical trajectories satisfying the equations of motion

Ṙ = ∂E 0 (R, P) ∂P , Ṗ = - ∂E 0 (R, P) ∂R ( 44
)
with the initial conditions

R(r, p; 0) = r , P(r, p; 0) = p , (45)
Let us now replace the quasiparticle-distribution function by N ν delta functions in phase space. Since the order of magnitude of ν 1 is dominated by -f ′ (E 0 ), it is clear that these delta functions should be distributed near the Fermi surface. To be more specific, we choose N ν points r i , p i in phase space which are distributed according to a probability density which is proportional to -f ′ (E 0 ), in such a way that for arbitrary but sufficiently smooth phase-space functions χ(r, p) the integral of χ(r, p) times the function f ′ [E 0 (r, p)] can be approximated by

d 3 rd 3 p (2π ) 3 χ(r, p)f ′ [E 0 (r, p)] ≈ -C Nν i=1 χ(r i , p i ) . (47)
Note that, if r i , p i are distributed in such a way, the same is true for R i (t) = R(r i , p i ; t), P i (t) = P(r i , p i ; t), since the quasiparticle energy

E i = E 0 [R i (t), P i (t)
] is a constant of the motion. In particular, defining y i (t) = y[R i (t), P i (t); t] and using Eq. ( 47), we can approximate the integral of an arbitrary function χ times the distribution function ν 1 as

d 3 rd 3 p (2π ) 3 χ(r, p)ν 1 (r, p; t) ≈ C Nν i=1 y i (t)χ[R i (t), P i (t)] .
(48) In other words, we have replaced ν 1 by

ν 1 (r, p; t) → C Nν i=1 y i (t)δ[r -R i (t)]δ[p -P i (t)] . (49)
According to Eq. ( 46), the equation of motion of the coefficients y i is reduced to

ẏi (t) = F [R i (t), P i (t); t] . ( 50
)
Above we assumed the function χ(r, p) to be sufficiently smooth. Of course, this causes some trouble if we want to calculate local quantities like the density or the current. For instance, we obtain

ρ 1ν (r) = C Nν i=1 y i (t) ξ i (t) E i δ[r -R i (t)] , (51)
where

ξ i (t) = h 0 [R i (t), P i (t)]
. This result makes sense only after the delta functions have been averaged over a volume containing a sufficiently large number of test particles in order to have a reasonable statistics. Supposing that this can be done, and supposing that V ext (r; t) and the phase φ 1 (r; t) are known, we can use the result for ρ 1ν in the explicit expression for F in order to obtain a system of N ν coupled first-order differential equations of the form (50) for the coefficients y i . This represents a tremendous simplification with respect to the original partial differential equation (38) in seven dimensions (r, p, and t). However, the phase φ 1 (r, t) is not known, but it has to be determined from the continuity equation (26). This is, again, very difficult. Hence, instead of solving the continuity equation exactly, we make an ansatz for φ 1 and determine the parameters by minimizing the violation of the continuity equation,

d 3 r( ρ1 + ∇ • j 1 ) 2 = min , (52)
the explicit expression for ρ1 + ∇ • j 1 being given by Eq. ( 29). The idea is to expand φ 1 on an appropriately chosen set of orthogonal functions ψ n ,

φ 1 (r; t) = N φ n=1 x n (t)ψ n (r) . (53)
Inserting this ansatz into Eq. ( 29), we see that the integral in Eq. ( 52) depends on x n and ẍn . At a given time, we regard x n and ẋn as given (e.g., at the moment when the perturbation is switched on, we know that x n = ẋn = 0). Hence, in order to have a minimal violation of the continuity equation, we have to minimize Eq. ( 52) with respect to ẍn by demanding

d dẍ n d 3 r( ρ1 + ∇ • j 1 ) 2 = 0 . ( 54
)
At this stage it turns out to be convenient to choose the basis functions ψ n such that they satisfy the orthogonality relation

d 3 r A 1 + gA 2 ψ n (r)ψ m (r) = δ nm . (55)
Then we obtain the following differential equation for the coefficients x n :

ẍn (t) = N φ m=1 a nm x m (t) + d 3 r A 2 ψ n (1 + gA) 2 V1ext -g∇ • j 1ν - ∆ 0 A ∇ • d 3 p (2π ) 3 ∆ 0 E 2 0 p m ν 1 , ( 56
)
where a is a time-independent matrix,

a nm = 2 m d 3 r A 2 ψ n (1 + gA) 2 2π 2 3 mp F + g ∇ • (ρ 0 ∇ψ m ) .
(57) Using Eq. ( 48) and integrating by parts, we can rewrite Eq. ( 56) in a more convenient form as

ẍn (t) = N φ m=1 a nm x m (t) + Nν i=1 b ni (t)y i (t) + vn (t) , ( 58
)
where b(t) denotes the matrix

b ni (t) = C m P i (t) • ∇ gA 2 ψ n (1 + gA) 2 + ∆ 0 E 2 i ∇ A∆ 0 ψ n (1 + gA) 2 R i (t) (59)
and the vector v is defined by

v n = d 3 r A 2 ψ n V 1ext (1 + gA) 2 . ( 60
)
Mainly for formal purposes, we note that also the equation (50) for the coefficients y i can be rewritten in matrix notation as

ẏi (t) = N φ n=1 [c in (t) ẋn (t) + d in (t)x n (t)] + f i (t) + Nν j=1 g ij (t)y j (t) , (61)
where

c in (t) = m P i (t) • ∇ ψ n 1 + gA - ∆ 0 E 2 i ∇ ∆ 0 ψ n 1 + gA R i (t) , (62)
d in (t) = m ξ i (t) E i (P i (t) • ∇) 2 ψ n m -(∇V 0 ) • ∇ψ n R i (t) , (63)
f i (t) = - P i (t) m • ∇ V 1ext 1 + gA - ∆ 0 E 2 i ∇ ∆ 0 V 1ext 1 + gA R i (t) , (64)
and

g ij (t) = -gC ξ j (t) E j P i (t) m • ∇ δ[r -R j (t)] 1 + gA - ∆ 0 E 2 i ∇ ∆ 0 δ[r -R j (t)] 1 + gA Ri(t)
.

In the latter equation, δ denotes a kind of "smeared" delta function which accounts for the averaging mentioned below Eq. (51). However, as mentioned above, Eqs. ( 61) and (65) will be used for formal purposes only. In practice, it is much faster to calculate ρ 1ν (r) after each time step on a discrete mesh, and to interpolate the stored values when performing the next time step for the coefficients y i . For the calculation of ρ 1ν (r), we replace the delta function in Eq. ( 51) by a Gaussian having a width d ρ .

In summary, the coupled system of partial differential equations, namely the transport equation for the distribution function ν 1 and the continuity equation for the phase φ 1 [Eqs. (38) and ( 26)], has been replaced by a coupled system of ordinary linear differential equations for the coefficients y i and x n [Eqs. (61) and (58)], which can formally be written as

d dt x(t) ẋ(t) y(t) = 0 1 0 a 0 b(t) d(t) c(t) g(t) x(t) ẋ(t) y(t) + 0 v(t) f (t)
. (66)

B. Trajectories of the test particles

In practice, the solution of the classical equations of motion for the test particles, Eqs. (44), faces us with some unusual features which are not present with the usual Newtonian equations of motion. Note that we are not dealing with ordinary particles but with Bogoliubov quasiparticles, which have some surprising properties. For instance, E i being a constant of the motion and

E 2 i = ξ 2 i + ∆ 2 0 (r i )
, it is evident that the energy ξ i cannot be conserved if the gap ∆ 0 depends on r. In particular, when a test particle with quasiparticle energy E i reaches the surface where ∆ 0 (r) = E i , it is reflected (Andreev reflection). During this reflection, the momentum P i stays almost constant, but the energy ξ i changes its sign (i.e., a particle is transformed into a hole or vice versa), such that the velocity v i = ∂E i /∂P i = (ξ i /E i )P i /m is reversed. As a consequence, the quasiparticle is reflected into the direction where it came from, which is very surprising if the incident angle is different from 90 • .

In order to find the test-particle trajectories numerically, it does not seem very efficient to start directly from Eqs. (44), since a small numerical error in the momentum of the order of δP/P ∼ ∆/ǫ F would immediately lead to a completely wrong behavior. It is therefore advantageous to make use of the variable ξ i , whose equation of motion reads

ξi = - ∆ 0 (R i ) E i P i m • ∇∆ 0 (R i ) . (67)
Solving this equation together with the equations for R i and P i , we can correct P i after each time step according to

P corr . i = P i |P i | 2mξ i + p 2 F (R i ) . (68)
In practice, the variable ξ i also allows us to introduce a very reliable method for determining the step size. Let us denote by ξ ′ i the result we obtain after one time step of size δt, and by ξ ′′ i the result we obtain after two time steps of size δt/2 each. Then the quantity δt|ξ ′ -ξ ′′ | is a measure for the numerical error and can be used for adapting the step size δt to the situation. It turns out that the step size has to become very small only during Andreev reflection. Now let us give some examples for typical test-particle trajectories. For that purpose, let us restrict ourselves to the most simple case which is a spherical harmonic trap,

V 0ext (r) = 1 2 mΩ 2 r 2 . ( 69
)
This potential defines the so-called trap units, i.e., energies are measured in units of Ω, temperatures in units of Ω/k B , lengths in units of l ho = /(mΩ), etc. In this example, due to spherical symmetry, not only the quasiparticle energy E, but also the angular momentum L = r × p of a test particle is a constant of the motion.

Within the local-density approximation (LDA) [START_REF] Houbiers | [END_REF]17], the density ρ 0 (r) has its maximum at the center of the trap and vanishes approximately (except for very small temperature effects) at the Thomas-Fermi radius R T F = 2µ/(mΩ 2 ). The gap ∆ 0 (r) has its maximum at the center of the trap, too, and goes to zero at some critical radius R c which is temperature dependent and determined by the equation T = T c (R c ). In order to avoid numerical problems arising from the infinite derivative of ∆ 0 (r) at r = R c , we convolute the LDA result for ∆ 0 (r) with a Gaussian of width d ∆ . In fact, this is more realistic than the LDA result since the exact solution of the BdG equations also leads to a gap ∆ 0 (r) which has an exponential tail [17,18,19]. As parameters we choose µ = 32 Ω, g = -2 l ho /m, and T = 1.4 Ω/k B . The corresponding number of atoms in the trap is approximately 1.7 × 10 4 . For these parameters quantum mechanical (BdG, QRPA) results are available for comparison. The width d ∆ is chosen such as to optimize the agreement with the BdG gap, which for the present parameters is achieved with d ∆ = l ho .

In Fig. 1 we show the corresponding gap ∆ 0 (r) as a function of the distance r from the center of the trap. From this figure it is evident that due to the condition E ≥ ∆ 0 (r), the relevant quasiparticles (having E < ∼ k B T = 1.4 Ω) are excluded from the region r < ∼ 4l ho . In addition to the gap, we display the potential V 0 (r) -µ, since the motion of a quasiparticle with given energy E and angular momentum L is also limited by the condition E 2 -∆ 2 0 (r) -L 2 /2mr 2 ≥ V 0 (r) -µ. It has been shown that also within the fully quantum-mechanical BdG theory the lowest-lying quasiparticle states are localized in this region [20]. In our example, the motion of the relevant quasiparticles is restricted to the region 4 < ∼ r/l ho < ∼ 8. Most of these quasiparticles will undergo Andreev reflection. Their trajectories are approximately described by an ellipse which is cut at the points where ∆ 0 (r) = E. If E ≪ ǫ F , the quasiparticle will move hence and forth on the same partial ellipse. Such trajectories with E = 0.1 Ω and 0.4 Ω are shown in the left panel of Fig. 2. However, if the quasiparticle energy is higher, the change in energy from ξ ≈ E to ξ ≈ -E (or vice versa) during the Andreev reflection results in a change of momentum which is no more negligible. Then, due to angular momentum conservation, the angle of reflection is slightly different from the angle of incidence, and the whole trajectory is precessing. An example for such a trajectory with E = 0.7 Ω is also shown in the left panel of Fig. 2. A completely different picture arises if the initial conditions are such that the quasiparticle does never reach the point where ∆(r) = E. Then the trajectory is just a precessing, slightly deformed ellipse, as shown in the right panel of Fig. 2 for the case of a trajectory with E = Ω. There is a striking analogy between these trajectories and the "glancing" orbits discussed, e.g., in Ref. [21] in the context of a superconducting cylinder which is coated by a normal-metal layer.

C. Distribution of test particles in phase space

In Sec. III A we supposed that one can generate a distribution of points r i , p i in phase space such that Eq. ( 47) is approximately satisfied for sufficiently smooth functions χ(r, p). In practice, this distribution is obtained in two steps. First we generate the coordinates r i , and in a second step the momenta p i .

The mean density of test particles at a certain point r is given by

n(r) = Nν i=1 δ(r i -r) , (70)
where δ denotes a smeared delta function in order to account for the averaging. Using Eq. ( 47), we conclude

n(r) = - 1 C d 3 p (2π ) 3 f ′ [E 0 (r, p)] ≡ w(r) C . ( 71
)
The algorithm for the generation of the coordinates r i is now very simple. First we look for the maximum w max of the function w(r). Defining P (r) = w(r)/w max , we obtain a function whose values lie between 0 and 1. Then we generate uniformly distributed random points r k in a volume which contains the whole system, and retain each point with the probability P (r k ), until the desired number of points, N ν , is reached. The formula (71) for the test-particle density n(r) can also be used for the determination of the normalization constant C. Integrating n(r) over space, we must recover the total number of test particles. This implies

C = 1 N ν d 3 r w(r) . ( 72
)
Now we turn to the distribution of the momenta p i . It is evident that the angular distribution of the momenta is isotropic, i.e., the interesting part of the problem is the distribution of the absolute values, p i = |p i |, which is, of course, directly related to the distribution of the energies ξ i . Let us define the mean number of test particles per energy and volume

n(r, ξ) = Nν i=1 δ(r i -r) δ(ξ i -ξ) . (73)
Again, with the help of Eq. ( 47), this becomes

n(r, ξ) = - 1 C mp ξ 2π 2 3 f ′ (E ξ ) , (74)
with p ξ = 2mξ + p 2 F (r), i.e., for given spatial coordinates r, the probability density for finding a particle at energy ξ is proportional to -p ξ f ′ (E ξ ). Such a distribution can be generated in the following way. Starting from random numbers z k which are uniformly distributed in the interval (0, 1), it is straight-forward to show that the energies

ξ k = T ln z k 1 -z k (75)
are distributed according to the probability density -f ′ (ξ). It is evident that negative energies with ξ < -ǫ F (r) have to be removed. Furthermore, it is preferable to cut the distribution at energies which lie too far away from the Fermi surface, e.g., |ξ| > 15T (the probability that this happens is less than 10 -6 ). The momenta p ξ are thus limited by p max = 30mT + p 2 F (r), and the function defined by P (ξ) = p ξ f ′ (E ξ )/p max f ′ (ξ) cannot become greater than 1 and can serve as a probability. If we retain each energy ξ k generated according to Eq. ( 75) with the probability P (ξ k ), the remaining energies are distributed according to the desired distribution. In order to give an illustration for the resulting distribution of test particles, we show in Fig. 3 the radial distribution of N ν = 10 5 test particles in a trap with the same parameters as in Fig. 1. In agreement with what we discussed in the preceding subsection, we see that the test particles are mainly located in the region 4 < ∼ r/l ho < ∼ 8, corresponding to the region where the system is mainly normal fluid. Due to the angular average the statistical fluctuations around the ideal distribution, Eq. ( 71), which is represented by the dotted line, are very small. We verified that, apart from the statistical fluctuations, our test-particle distribution stays constant, which is a good numerical test of both the initial test-particle distribution and of the test-particle trajectories.

D. Initial condition

In the linear response regime, as the name implies, the response to a time-dependent perturbation of the form V 1 (r; t) = V1 (r)f (t), with an arbitrary time dependence f (t), can be obtained as convolution of f (t) with the response to a delta function in time. It is therefore sufficient to study perturbations of the form

V 1ext (r; t) = V1 (r)δ(t) . (76)
We thus set the inhomogeneous terms in Eq. (66) to v(t) = v δ(t) and ḟ (t) = f δ(t), respectively, v and f being defined analogously to Eqs. ( 60) and (64) but with V 1ext replaced by V1 . Assuming that the system was in equilibrium before this perturbation, we may ask the question: What are the values of the coefficients y i and x n immediately after the perturbation, i.e., at infinitesimally small t > 0? This question can be answered exactly, since during the infinitesimal period where the perturbation is active, the matrix in Eq. (66) can be regarded as time-independent. Integrating Eq. (66) over time from -t 0 to t 0 , we obtain in the limit

t 0 → 0 lim t0→0 x(t 0 ) ẋ(t 0 ) y(t 0 ) = v 0 cv + f . (77)
Let us now assume that the function V1 lies in the space spanned by the functions ψ n . Then it is evident that the corresponding linear combination is given by the coefficients vn , i.e.,

V1 (r) = N φ n=1 vn ψ n (r) . (78)
Note that the functions ψ n do not necessarily have to have this property. For example, we could define a basis of functions satisfying the orthogonality relation ( 55) and maybe even a suitably defined completeness relation if N φ → ∞, but which all vanish identically outside the superfluid region, i.e., in the region where ∆ 0 = 0 (and A = 0). Eq. ( 78) would then be satisfied inside the superfluid region, but not outside. Hence, it is an additional requirement for the choice of the functions ψ n . Combining Eqs. ( 77) and (78), we find

lim t0→0 φ 1 (r; t 0 ) = 1 V1 (r) . (79)
Eq. ( 78) also leads to a simplification of the initial value of the coefficients y i and the quasiparticle-distribution function. Using the explicit expressions for the matrix c and the vector f [Eqs. ( 62) and ( 64) with V 1ext replaced by V1 ], we obtain from the third line of Eq. ( 77)

lim t0→0 y i (t 0 ) = N φ n=1 c in vn + fi = p i m • ∇ N φ n=1 vn ψ n -V1 1 + gA - ∆ 0 E 2 i ∇ ∆ 0 ( N φ n=1 vn ψ n -V1 ) 1 + gA ri . (80)
As a consequence, if Eq. ( 78) is satisfied, the initial values of the coefficients y i vanish, which implies

lim t0→0 ν 1 (r, p; t 0 ) = 0 . (81)
In fact, the simple result of this subsection, which is summarized in Eqs. ( 79) and (81), could have been anticipated without any calculation. The effect of a perturbation of the form (76) is to give a particle at position r a kick

δp = -dt∇V 1ext (r; t) = -∇ V1 (r) . (82)
Since this kick does not depend on the momentum of the particle, the local Fermi sphere is shifted as a whole, there is no change in density and no Fermi surface deformation. Within the present theoretical framework, Cooper pairs are not broken either, they just acquire a center of mass momentum. Thus, the distribution function in the local rest frame stays unchanged (ν 1 = 0), and the collective velocity is given by

v coll = -( /m)∇φ 1 = -(1/m)∇ V1 .
Note, however, that in reality a perturbation which has the form of a short pulse would lead to much more complicated effects (e.g., pair breaking). Since our semiclassical description requires that the time dependence of the perturbation is slow, our formal result for a deltalike excitation becomes physically meaningful only after it has been convoluted with a function f (t) which varies slowly in time. In other words, we can only calculate the low-frequency part of the response function.

IV. FIRST RESULTS

In this section we will discuss first numerical results which have been obtained using the test-particle method. Our intention here is to see whether this method is in principle capable to describe the most important features of collective excitations in superfluid trapped Fermi gases. To that end, we will study the quadrupole excitation of a spherical system, which is excited by

V1 (r) = α mΩ(2r 2 z -r 2 x -r 2 y ) (83)
(the factor mΩ has been introduced in order to make the coefficient α dimensionless).

For practical purposes, we will make an additional approximation: We will restrict our ansatz for the phase, Eq. ( 53), to only one or two functions ψ n . It is clear from rotational symmetry that in the case of a quadrupole excitation of the form (83) the most general form the phase can have is

φ 1 (r) = Φ(r)[2r 2 z -r 2 x -r 2 y ] , (84)
such that that the functions ψ n can be written as

ψ n (r) = Ψ n (r)[2r 2 z -r 2 x -r 2 y ] . (85)
It is known from superfluid hydrodynamics that at zero temperature the velocity field is essentially linear in the coordinates, i.e., the function Φ(r) is almost constant. As a first guess we will assume that this is still true at nonzero temperature, and hence we will take only one single function (N φ = 1) in the ansatz (53) for the phase, Ψ 1 = const. The proportionality constant will be determined from the normalization condition (55). Such a restricted ansatz means of course that the continuity equation will not be exactly satisfied in the superfluid region (remember that outside the superfluid region the phase has no effect whatsoever). We will therefore improve this initial ansatz by including a second function (N φ = 2) which allows to modulate Φ(r) in the superfluid region.

The first idea one might have is to use for Ψ n (r) polynomials in r 2 and to orthogonalize the resulting functions ψ n . However, it turns out that this leads to numerical instabilities due to the fast growing of the resulting polynomials outside the superfluid region. Let us explain this effect in some more detail. As seen from the transport equation for the quasiparticle-distribution function, the phase φ 1 outside the superfluid region enters directly the dynamics of ν 1 . Although the net effect of the phase and of the quasiparticle motion should be independent of the choice of φ 1 outside the superfluid region, each of these contributions depends on this choice. If φ 1 changes too rapidly, the numerical solution of the equation of motion for the coefficients y i becomes less accurate and the cancellation of the two effects does not work any more.

We therefore have to look for functions Ψ n which are linearly independent inside the superfluid region, but which do not grow outside. Here we will choose the functions Ψ1 (r) = 1 and Ψ2 (r) = [1 -ϕ(r)] 2 , where ϕ(r) is the function defined in Eq. (25). The function Ψ2 (r) has its maximum in the center of the trap and goes smoothly to zero at the boundary of the superfluid region. From Ψ1 and Ψ2 the functions Ψ 1 and Ψ 2 are determined according to the orthogonality condition (55) with the help of the Gram-Schmidt orthogonalization method. As we will see, the results obtained with N φ = 1 and N φ = 2 are very similar and we therefore claim that they would not change a lot if we included additional functions.

Let us now present the results. As in the examples shown in the preceding section, we consider a spherical harmonic trap with µ = 32 Ω, containing approximately 1.7 × 10 4 atoms. The resulting density profile ρ 0 (r) is shown in Fig. 4 as the dashed line. The critical temperature within LDA is T c = T c (r = 0) ≈ 3.9 Ω/k B . As before, the LDA result for the gap ∆ 0 (r) is convoluted with a Gaussian having a width d ∆ = l ho . We will study the quadrupole mode for three different temperatures, T /T c = 0.2, 0.4, and 0.6. The equilibrium gap ∆ 0 (r) for these three temperatures is also displayed in Fig. 4.

After the system is excited, its shape will oscillate. A measure for this quadrupole deformation is the ratio

2r 2 z -r 2 x -r 2 y r 2 0 , (86)
where r 2 0 denotes the mean square radius in equilibrium, which in the present case has the value r 2 0 ≈ 23 l 2 ho . In the linear response, the quadrupole deformation is of course proportional to the strength of the perturbation, and we therefore divide our results by this strength [denoted α in Eq. ( 83)]. In our simulation we use N ν = 10 5 test particles, the width of the Gaussians used for smearing ρ 1ν (see Sec. III A) is set to d ρ = l ho . In Fig. 5 we display the time dependence of the quadrupole deformation after the perturbation for the three temperatures mentioned before. The corresponding spectra, obtained by Fourier transformation, are shown in Fig. 6. The results for the two cases N φ = 1 and N φ = 2 are displayed as dashed and solid curves, respectively. In all cases the two curves are in reasonable agreement, such that we can say that the use of N φ = 2 independent functions in the ansatz for the phase is sufficient.

We see that the temperature dependence of the spectrum is highly non-trivial. At low temperatures (upper panel of Fig. 6), we see essentially the hydrodynamic quadrupole mode, which at zero temperature lies at ω = √ 2Ω [22,23,24] and which is now damped as a consequence of its coupling to the normal component. At higher temperatures (middle of Fig. 6), a second peak builds up in the spectrum, corresponding to the quadrupole mode in the normal phase, which lies slightly above ω = 2Ω [25] (for the present set of parameters, its frequency is ω ≈ 2.2 Ω [6]). As the temperature approaches T c (lower panel of Fig. 6), the strength contained in this second peak increases, while the hydrodynamic mode, whose frequency is slightly shifted downwards, disappears. These findings are in good agreement with quantum mechanical QRPA calculations [6].

We note that the damping width of the hydrodynamic mode at low temperature is now comparable with that found within the QRPA and much stronger than that found in our previous work [10], where we replaced the gap ∆ 0 (r) by a constant. The reason is in fact very simple: With a constant gap, the fraction ρ n /ρ 0 of the normal component is independent of r, whereas in the case of an r-dependent gap the normal component in the outer part of the system is already important at very low temperatures [3].

As the temperature T approaches T c , the quadrupole mode of the normal phase (that at ω = 2.2 Ω) becomes undamped, as it is the case within the QRPA. However, even though collisions are strongly suppressed at these low temperatures, it should be kept in mind that the collision term, which is neglected in the present work, is non-zero and its inclusion would lead to a finite lifetime of this oscillation, too.

Finally, let us compare our semiclassical results more quantitatively to quantum mechanical QRPA results. In Fig. 7 we show the QRPA result of Ref. [6] for the quadrupole excitation spectrum (dotted line) together with the semiclassical result we obtain with the same parameters (solid line). As one can see, the total normalization and the relative weights of the two peaks are in reasonable agreement. Also the widths of the QRPA peaks are well reproduced by the semiclassical calculation. The main differences are that within the semiclassical calculation the two peaks lie a bit too high and that they are not as well separated as within the QRPA. A comparison of the two curves for N φ = 1 and N φ = 2 in the T /T c = 0.4 case shown in the middle of Fig. 6, whose parameters are quite close to those of Fig. 7, suggests that the latter effect might be partly due to the restricted ansatz for the phase. However, as one can deduce from the irregular structure of the QRPA spectrum, even in a system with 32000 atoms shell effects, i.e., effects which depend on the discrete single-particle spectrum, are still quite pronounced. It is clear that such effects cannot be reproduced within a semiclassical calculation. In this sense the agreement between the two spectra is very satisfactory, in particular since one can assume that the shell effects decrease with increasing number of particles.

V. CONCLUSIONS

In this paper, we developed a numerical test-particle method for solving the semiclassical transport equations for an ultracold trapped Fermi gas in the BCS phase in the collisionless limit. These transport equations take into account the coupling between the dynamics of the Cooper pairs (superfluid component) and the thermally excited Bogoliubov quasiparticles (normal component). We developed the method for the case of small deviations from equilibrium, so that the test-particle trajectories can be calculated in the equilibrium state. Since the testparticles describe Bogoliubov quasiparticles rather than real particles, the trajectories have very unusual properties compared with the trajectories one has to deal with when applying the test-particle method to the normal Vlasov equation. Our test particles can have the character of particles as well as holes, depending on whether their energy ξ is positive or negative, and they can also be transformed from the one into the other if they hit the region where the gap ∆ becomes larger than their quasiparticle energy E (Andreev reflection). Another complication as compared with the normal Vlasov equation is that the dynamics of the quasiparticles is coupled to the collective motion of the superfluid component, which is described by the phase φ of the order parameter. This phase has to be determined simultaneously with the evolution of the quasiparticle-distribution function by solving the continuity equation. In the present work, we make an ansatz for φ with time-dependent coefficients, leading to an approximate solution of the continuity equation.

As a first application, we calculated the response of a gas trapped in a spherical trap to a delta-like perturbation of quadrupole form. After this perturbation, the shape of the gas shows a damped oscillation. At low temperatures, this oscillation is just the hydrodynamic quadrupole mode which is damped by its coupling to the normal component. With increasing temperatures, the extension of the normal component increases, and, as a consequence, the normal component can perform its own quadrupole oscillation. Since the frequency of the quadrupole mode in the normal collisionless Fermi gas is higher than that of the hydrodynamic mode, this leads to a two-peak structure in the response function. As the temperature approaches T c , the strength of the hydrodynamic mode disappears and only the normal mode survives.

The next step will be to apply the method presented here to more realistic cases, namely to the axial and radial breathing modes of a gas in a cigar-shaped trap containing a larger number of particles. In fact, the deformation and the large particle number do not pose a big problem, which is one of the main advantages of the present method as compared with quantum mechanical QRPA calculations. Another possible application of the method is to study the dynamics of a vortex, where already the equilibrium situation is characterized by a non-vanishing phase of the order parameter.

However, there are still a number of unsolved problems and possible improvements of the method. First of all, the collision term [START_REF] Vollhardt | The Superfluid Phases of Helium 3[END_REF] should be included, which is an additional source of damping of collective oscillations. As mentioned in the introduction, the possibility to include collisions is an important advantage of the present method as compared with the QRPA, where collision effects cannot be taken into account since this would necessitate to include four-quasiparticle excitations. Second, from a fundamental point of view, the fact that the continuity equation is only approximately fulfilled is of course unsatisfactory and one should think about another numerical method for solving the continuity equation. Finally, one might ask the question how the present theory can be extended to the strongly interacting regime. Unfortunately, this question is up to now completely open, since in this regime thermal fluctuations of the order parameter, which are not contained in the BdG equations, play a crucial role (see, e.g., Ref. [26]).

FIG. 1 :

 1 FIG.1: Gap ∆0(r) (solid line) and potential V0(r)-µ (dashed line) for the case of a spherical trap with frequency Ω, chemical potential µ = 32 Ω, coupling constant g = -2 l ho /m, and temperature kBT = 1.4 Ω. ∆0 and V0 -µ are in units of Ω, r is in units of the oscillator length l ho . Roughly speaking, these two curves determine the classically allowed region for a quasiparticle with given energy E.

1 FIG. 2 :

 12 FIG. 2: Four examples of quasiparticle trajectories in a trap with parameters given below Fig. 1. The three trajectories shown in the left panel belong to quasiparticles with E = 0.1 Ω, 0.4 Ω, and 0.7 Ω, respectively. The trajectory displayed in the right panel belongs to a quasiparticle with E = Ω.

FIG. 3 :

 3 FIG.3: Radial distribution of 10 5 test particles in a trap with parameters given below Fig.1, counted in bins of size δr = 0.1 l ho . For comparison, the dotted curve represents the ideal distribution according to Eq. (71).

FIG. 4 :

 4 FIG.4: Density profile ρ0(r) (dashed line) and gap ∆0(r) (solid lines) in a spherical harmonic trap containing 1.7 × 10 4 atoms (µ = 32 Ω, interaction strength g = -2 l ho /m). The gap is displayed for three different temperatures, T /Tc = 0.2, 0.4, and 0.6, while the density profile is practically independent of T .

FIG. 5 :

 5 FIG.5: Time dependence of the quadrupole deformation after a delta-like perturbation at t = 0. The parameters are the same as in Fig.4, the three panels correspond, from top to bottom, to T /Tc = 0.2, 0.4, and 0.6. The dashed and solid lines correspond to N φ = 1 and N φ = 2, respectively.

FIG. 6 :FIG. 7 :

 67 FIG. 6: Fourier transforms of the quadrupole responses shown in Fig. 5.

Acknowledgments

The author wishes to thank P. Schuck for numerous fruitful discussions and the critical reading of the manuscript.


Gorky Shaw
email: gorky.shaw@ulg.ac.be.

R B G Kramer

N M Dempsey

K Hasselbach

A Scanning Hall Probe Microscope for high resolution, large area, variable height Magnetic Field Imaging

We present a Scanning Hall Probe Microscope operating in ambient conditions. One of the unique features of this microscope is the use of the same stepper motors for both sample positioning as well as scanning, which makes it possible to have a large scan range (few mm) in x and y directions, with a scan resolution of 0.1 µm. Protocols have been implemented to enable scanning at different heights from the sample surface. The z range is 35 mm. Microstructured Hall probes of size 1-5 µm have been developed. A minimum probe-sample distance < 2 µm has been obtained by the combination of new Hall probes and probe-sample distance regulation using a tuning fork based force detection technique. The system is also capable of recording local B(z) profiles. We discuss the application of the microscope for the study of micro-magnet arrays being developed for applications in micro-systems.

I. INTRODUCTION

Scanning Hall Probe Microscopy (SHPM) is a popular local magnetic characterization technique [START_REF] Chang | Scanning Hall probe microscopy[END_REF][START_REF] Oral | Real-time scanning hall probe microscopy[END_REF][START_REF] Oral | Scanning Hall probe microscopy of superconductors and magnetic materials[END_REF][START_REF] Howells | Scanning Hall probe microscopy of ferromagnetic structures[END_REF][START_REF] Sandhu | Room temperature scanning Hall probe microscopy using GaAs/AlGaAs and bi micro-hall probes[END_REF]. Among common non-invasive magnetic imaging techniques, SHPM is particularly useful due to a number of advantages. In comparison to Magnetic Force Microscopy (MFM) and Magneto-Optical Imaging (MOI), SHPM enables direct quantitative mapping of magnetic field. MFM is more challenging for quantitative measurements as characterization and modeling of the magnetization of the tip are rather difficult [START_REF] Bending | Scanning Hall probe microscopy of vortex matter[END_REF]. In comparison to MOI, SHPM works in a much larger magnetic field range, and offers better spatial resolution in the case when the sample is not completely flat [START_REF] Kirtley | Fundamental studies of superconductors using scanning magnetic imaging[END_REF]. Though SHPM involves longer image acquisition time than MOI unlike in a scanning probe system, in MOI the whole image is acquired at one time, and involves acquisition times of a few seconds at most) and offers worse spatial resolution compared to MFM, the unique compromise between spatial resolution and field sensitivity makes it a very useful tool for convenient non-destructive quantitative characterization over wide temperature and field ranges. Over the years, efforts have been made to improve the spatial resolution and scan range [START_REF] Gregory | A scanning Hall probe microscope for large area magnetic imaging down to cryogenic temperatures[END_REF][START_REF] Dinner | Cryogenic scanning Hall-probe microscope with centimeter scan range and submicron resolution[END_REF][START_REF] Cambel | Large-scale high-resolution scanning hall probe microscope used for mgb 2 filament characterization[END_REF][START_REF] Tang | An interchangeable scanning hall probe/scanning squid microscope[END_REF], temperature range [START_REF] Özgür Karcı | An ultra-low temperature scanning Hall probe microscope for magnetic imaging below 40 mK[END_REF], and magnetic field resolution [START_REF] Shimizu | Scanning Hall probe microscopy with high resolution of magnetic field image[END_REF] of SHPM.

In this paper we report on a room temperature SHPM system developed to characterize arrays of micro-magnets. SHPM is the tool of choice for quantifying the stray fields produced by micro-magnets, which have many applications in bio-medical studies [START_REF] Zanini | Autonomous micro-magnet based systems for highly efficient magnetic separation[END_REF][START_REF] Brunet | Evolutionary conservation of early mesoderm specification by mechanotransduction in bilateria[END_REF][START_REF] Pivetal | Magnetic nanoparticle DNA labeling for individual bacterial cell detection and recovery[END_REF] and MEMS [START_REF] Dempsey | Hard magnetic materials for MEMS applications[END_REF]. The present SHPM system demonstrates the usefulness of a combination of concepts, tuning fork feedback based height control, monolithic electronics for regulation and Hall signal read out, scanning with stepper motors, compared to microscopes presented before [START_REF] Gregory | A scanning Hall probe microscope for large area magnetic imaging down to cryogenic temperatures[END_REF][START_REF] Dinner | Cryogenic scanning Hall-probe microscope with centimeter scan range and submicron resolution[END_REF][START_REF] Cambel | Large-scale high-resolution scanning hall probe microscope used for mgb 2 filament characterization[END_REF][START_REF] Tang | An interchangeable scanning hall probe/scanning squid microscope[END_REF][START_REF] Özgür Karcı | An ultra-low temperature scanning Hall probe microscope for magnetic imaging below 40 mK[END_REF][START_REF] Shimizu | Scanning Hall probe microscopy with high resolution of magnetic field image[END_REF][START_REF] Kustov | Magnetic characterization of micropatterned Nd-Fe-B hard magnetic films using scanning Hall probe microscopy[END_REF]. Microstructured Hall probes with active area of size 1-5 µm have been fabricated for use in the microscope. The characteristics of the microscope include: large in-plane scan range (up to a few mm, limited by sample topography) with fine scan resolution (step resolution of the motors is 0.1 µm, which is, in principle, the minimum possible scan resolution. However, the spatial resolution of the magnetic field distribution obtained in this technique is limited by the size of the Hall probes as the detected magnetic field is the average field over the active area of the probe) and large z range, allowing sample-probe distances from < 2 µm to 35 mm, with a step resolution of 0.2 µm, as well as large magnetic field detection (fields up to ∼ 1 T have been measured) with high field resolution (100 µT).

Figures 1(a) and 1(b) show a schematic and an optical image, respectively, of the microscope. The microscope is mounted on a vibration-isolation table. Vibrations from the floor are filtered by a concrete block supported by rubber dampers filled with compressed air flowing viscously between the dampers and air reservoirs. The microscope is mounted on a copper block that rests on the concrete block. The microscope consists of two separate parts, viz., a mobile part which consists of a sample holder placed on an XYZ stage, and a fixed part which consists of the probe mounted on a tuning fork. This second part includes piezoelectric elements for z-regulation, tuning fork excitation, and all electrical connections involved in the measurement. Note that it is important that both of these parts rest on the same base (in this case the copper block) as otherwise additional and uncontrolled tilt between the sample and probe planes would arise which would lead to unnecessary complications in sample-probe distance regulation. The microscope is covered with a 50 × 50 × 60 cm 3 opaque enclosure to minimize the influence of ambient light. Figure 2 shows a schematic of the complete setup. A fast lockin-amplifier (Zurich Instruments HF2LI, henceforth referred to as 'the HF2LI') [19] serves as the primary electronics component of the microscope. It is used to provide the excitation signal as well as to measure the current from the tuning fork. It is also used to provide the Hall probe excitation current as well as to measure the generated Hall voltage. Finally, it also provides the PID control for the z-piezo stack mentioned (explained below). All control and measurement operations are performed through custom-built reliable and user-friendly LabVIEW programs developed specifically for the setup. in X and Y directions is 50 mm and the minimum step size is 0.1 µm. This minimum step size is precise enough for the primary purposes the microscope is built for. In our microscope, the stepper motors themselves perform the scanning, additional piezo scanners are superfluous. As a result, in principle, the maximum achievable scan range is the same as the travel range of the X/Y motors = 50 mm. However, in practice, the effective scan range is limited by the properties of the sample being measured. Most often, irregularities of the sample surface mean that the probe-sample distance regulation (discussed in Section II A 3) can be maintained only over a smaller distance, thus limiting the scan range. We have been able to perform scans over areas as large as 2.5 mm × 2.

Probe stage

The fixed probe stage consists of the Hall probe glued to a commercial quartz tuning fork [23] which is in turn glued to a probe plate (Cu coated fibre glass plate). Electrical connections to the Hall probe and tuning fork are anchored on the probe plate. The tuning fork is mechanically excited through a 'thickness-mode' piezoceramic element below the tuning fork (labelled 'Excitation piezo' in Fig. 1(a)) [START_REF] Veauvy | Scanning -superconduction quantum interference device force microscope[END_REF][START_REF] Hykel | MicroSQUID Force Microscopy in a Dilution Refrigerator[END_REF]. The excitation signal is provided from the first output of the HF2LI (cf. Figure 2). This assembly is attached to a 'z-piezo stack' extension piezoelectric element [26] which is used for regulation of the probesample distance. The z-piezo stack has a vertical extension range ∼ 32 µm (with input voltage range 0-120 V). Voltage on the z-piezo stack is applied through one of the auxiliary outputs of the HF2LI via a low noise voltage amplifier with 10x amplification (cf. Figure 2). The piezo stack is mounted on a base block that can be tilted in order to adjust the angle between the Hall probe and the sample surface.

Probe-sample distance regulation

One of the primary challenges in a scanning probe setup is the ability to approach the probe as close as possible to the sample surface without crashing the probe into the sample. This is achieved by real-time control of the probe-sample distance by regulating the extension of the z-piezo stack during a scan using the amplitude of vibration of the tuning fork for feedback. Tuning fork feedback based probe-sample distance control was first proposed by Khaled Karrai and Robert D. Grober [START_REF] Karrai | Piezoelectric tip-sample distance control for near field optical microscopes[END_REF][START_REF] Karraï | Piezo-electric tuning fork tip-sample distance control for near field optical microscopes[END_REF] and has since been widely used in near-field scanning optical microscopy (NSOM). This technique has been implemented in Scanning SQUID [START_REF] Veauvy | Scanning -superconduction quantum interference device force microscope[END_REF][START_REF] Hykel | MicroSQUID Force Microscopy in a Dilution Refrigerator[END_REF] and Scanning Hall probe [START_REF] Dede | Scanning Hall Probe Microscopy (SHPM) Using Quartz Crystal AFM Feedback[END_REF] Microscopes. The voltage on the z-piezo stack is regulated by a fast closed-loop Proportional-Integral-Derivative (PID) control in the HF2LI PID settings have been optimized based primarily on the sample surface type. Typical values for samples with hard surfaces (e.g., thermomagnetically patterned micromagnets) are: P = -0.05, I = -50 s -1 , D = +1 µs. For samples with soft surfaces (e.g., topographically patterned micromagnets embedded in polymer) typical values of the parameters are: P = -0.5, I = -500 s -1 , D = +4 µs. These settings correspond to an integral time constant (P/I) = 1 ms. The PID output is updated at 100 kHz. The input of the PID is the amplitude of vibration of the tuning fork (cf. Figure 2). Current proportional to the amplitude of vibration is generated in the tuning fork, which is amplified with a current/voltage converter of gain 10 7 V/A, and passed through a high-pass filter (cut-off frequency 8 kHz). The amplified voltage is demodulated by the HF2LI and the phase signal is fed into a phase locked loop (PLL). The PLL ensures that the tuning fork oscillation frequency always remains in the vicinity of the resonance frequency. The excitation frequency is adjusted in order to maintain the tuning fork at resonance (cf. Figure 2). The resonance frequency f R of the free tuning fork is ∼ 32.768 kHz. However, it depends strongly on environmental conditions. Figure 3(a) shows a typical frequency response of a tuning fork glued to the probe plate.

Note that f R is reduced to ∼ 31.45 kHz. Figure 3(b) shows the frequency response after a probe is glued to the tuning fork. In this case f R is further reduced to ∼ 22.76 kHz.

f R for a tuning fork with a probe attached varies in the range 14-23 kHz. Note that the amplitude of oscillation of the tuning fork also reduces with increase in load. This is not obvious from Figures 3(a) and 3(b) as the excitation voltage in each case is adjusted such that the amplitude at resonance is always ∼ 1 V.

Further, the amplitude of oscillation monotonically decreases on approaching the sample surface. This is the basis of the probe-sample distance regulation [START_REF] Karrai | Piezoelectric tip-sample distance control for near field optical microscopes[END_REF][START_REF] Karraï | Piezo-electric tuning fork tip-sample distance control for near field optical microscopes[END_REF]. Figure 3(c) shows a typical approach curve for a loaded tuning fork (i.e., after a probe is mounted on the tuning fork glued to a probe plate). The x-axis shows the distance of the tuning fork from the sample surface (sample surface at z = 0). Note from Figure 3(c) that as the probe approaches the sample surface, the resonance frequency increases before decreasing, while the amplitude of oscillation decreases monotonically. Hence it is possible to maintain a fixed desired probe-sample distance if the amplitude of oscillation can be regulated. This is implemented via the PID control of the extension of the z-piezo stack using the amplitude of oscillation of the tuning fork as input as discussed above. In general, 70% of amplitude of oscillation of the tuning fork when it is far away from the sample surface is used as the desired set-point for probe-sample contact. For samples with rough and/or soft surfaces, 75-80% is set. From Figure 3(c), z ≈ 470 nm for 70% and 500 nm for 80%. This is the probe-sample distance which is maintained during a scan in the regulation mode (cf. Section II C 2).

B. Magnetic field detection

Hall probes

The Hall probes consist of heterostructures of GaAs/AlGaAs two-dimensional electron gas (2DEG) materials. The two dimensional electron gas is about 50 nm below the surface.

The mesa is first etched with a chemical etch by 100 nm in order to structure the Hall cross.

Then a multilayer of Ni/Au/Ge is evaporated with an e-gun evaporator (Plassys). Rapid thermal annealing allows the diffusion of the AuGe to the two-dimensional electron gas, thus forming the Ohmic contacts with the 2DEG. An additional deep etch ∼ 10 µm is made in order to place the Hall cross exactly at the corner of the wafer. Finally, the Hall probes are cut with a diamond dicing machine. The electron sheet density n is ∼ 3 × 10 12 cm -2 (calculated from the relation n = 1/(ke), where e is the electronic charge and k is the Hall sensitivity = 200 Ω/T, cf. Section II B 2 below). The mobility of the electron gas measured at 4.2 K is ∼ 30000 cm 2 /(Vs). Figure 4 Though Hall probes of smaller size have been prepared with GaAs/AlGaAs, smaller sizes are also associated with lower field sensitivity [START_REF] Chang | Scanning Hall probe microscopy[END_REF][START_REF] Sandhu | Room temperature scanning Hall probe microscopy using GaAs/AlGaAs and bi micro-hall probes[END_REF][START_REF] Hicks | Noise characteristics of 100 nm scale GaAsal[sub x]ga[sub 1-x]as scanning hall probes[END_REF], hence the size of the present probes is suitable for our purposes.

Probe calibration and magnetic field measurement

The Hall probe is excited with an ac signal from the second output of the HF2LI through a 100kΩ resistor in series. The Hall excitation current is set in the range 10-50 µA (pk-pk).

The frequency of excitation is set at 522 Hz. This was chosen empirically to correspond to minimum noise in the output signal. The amplified and filtered Hall voltage is detected at the second input of the HF2LI (cf. Figure 2). A number of Hall probes have been used in the measurements. All the Hall probes show remarkably consistent magnetic field sensitivity.

The Hall probes were calibrated using a pre-calibrated copper coil magnet (16.83 G/A, or 1.683 mT/A) [31]. For calibration, a Hall probe is placed at the center of the magnet and the Hall voltage is measured for a set of values of current flowing through the solenoid (corresponding to magnetic field values in the range -8 mT to +8 mT, in the magnet). Where I 24 is the current applied across the sides 1 and 3 of the Hall cross and V 13 is the Hall voltage (in-phase component of the demodulated signal) measured across the sides 2 and 4.

From Figure 4(d) it is clear that the Hall resistance varies very linearly with magnetic field.

From the linear fit, the Hall sensitivity (k) = is found to be 200 ± 2 Ω/T for all the used probes, which is consistent with the fact that Hall probe sensitivity is nearly independent of the Hall cross area [START_REF] Kirtley | Fundamental studies of superconductors using scanning magnetic imaging[END_REF]. In our Hall probes the doping level is rather high which results

in the relatively lower Hall sensitivity compared to some other reported values at room temperature (300 K) for GaAs/AlGaAs heterostructures (3000 Ω/T in Ref. ) however, it is sufficient for our purposes. We believe this also ensures that conduction always takes place preferentially in the 2DEG as there are many charge carriers available. Indeed, the Hall sensitivity was measured at room temperature as well as down to 50 mK and is found to be nearly constant in this temperature range (variation smaller than 5%).

Once the sensitivity of the Hall probe is known, the measured Hall voltage V 13 can be simply converted to magnetic field B z as,

B z = V 13 /I 24 /k (2)
Or,

B z = R Hall /k (3)
Note that there is always some offset voltage, i.e., a non-zero measured Hall voltage even in zero field. This can appear due to various reasons, including misalignment of the voltage leads and/or leakage between the leads. The offset voltage is measured before each measurement and the magnetic field value is determined by subtracting this offset signal from the measured Hall signal.

The magnetic field resolution is limited by the noise in the output signal. The noise depends on the active area of the probe and is larger for a smaller probe. In order to obtain understanding of the sensitivity limitations of our microscope let us compare the observed noise level to the Johnson-Nyquist thermal noise which limits the sensitivity of measurements [START_REF] Oral | Real-time scanning hall probe microscopy[END_REF][START_REF] Geim | Ballistic hall micromagnetometry[END_REF][START_REF] Novoselov | Submicron probes for hall magnetometry over the extended temperature range from helium to room temperature[END_REF]:

V johnson = 4k B T R (4)
Where k B is the Boltzmann constant, T is the temperature and R is the resistance of the Hall sensor. At room temperature (T = 300 K), using R = 120 kΩ (maximum contact resistance for the smallest probe, including the series resistance 100 kΩ), (± 0.5 • , measured optically), which is small enough to assume that the measured signal the component of the magnetic field is perpendicular to the sample surface [31]. Note that α ∼ 5 • is usually good enough to ensure that the micro-bonded wires do not touch the sample surface. However, a larger angle is needed for safety of the probe itself. Since the Hall cross is rather close to the edge of the Hall bar, scanning over large areas of samples with large topographic variations leads to significant bruising of the probe. While scanning with a smaller angle α, in a number of probes, the contact closest to the edge got disconnected due to scratches induced during scanning, rendering the probe unusable. The angle β is kept as small as possible (∼ 1 • ). The distance h of the Hall cross from the sample surface is then given by,

V
h = L 1 sinα + L 2 sinβcosα (5)
Using L 1 = L 2 = 8 µm (15 µm), α = 9 • and β = 1 • in Equation 5, we get, h ≈ 1.4 µm (2.6 µm) for the smallest (largest) probe, which is the minimum possible probe-sample distance that can be achieved without endangering the probe.

The distances of the Hall cross from the Hall bar edge are comparable to some other reported SHPM systems (13 µm [START_REF] Oral | Real-time scanning hall probe microscopy[END_REF], 20 µm [START_REF] Tang | An interchangeable scanning hall probe/scanning squid microscope[END_REF]) and are significantly smaller than that in the probes used in the previous setup [START_REF] Kustov | Magnetic characterization of micropatterned Nd-Fe-B hard magnetic films using scanning Hall probe microscopy[END_REF] (100 µm). As evident from the discussion above, the smaller distance in the present Hall probes has allowed achieving a much smaller probesample distance in the present setup ∼ 2 µm as compared to ∼ 25 µm in the previous setup).

Scan protocol: regulation and flyover modes

The scanning plane is designated as the (x, y) plane and hence the measured magnetic field component is along the z-axis (B z ). As discussed in Section II A 3, the probe is mounted on a quartz tuning fork allowing topographic feedback to maintain constant probe -sample distance during a scan. After setting the desired starting (x, y) position, the sample is approached towards the probe with active topographic feedback via the PID control, till the desired probe-sample distance is reached (< 500 nm, determined by the set amplitude of oscillation of the tuning fork, usually 70% of the amplitude of oscillation when the tuning fork is far away from the sample surface, cf. Section II A 3). The fast scan direction is along the y-axis, i.e., the long side of the Hall bar (cf. Figure 4). Scans are performed in the y-direction after setting the desired x before each scan in y. While making a step in y during the scan, the probe -sample distance is continuously regulated through the fast PID control. At each y-step, the magnetic field signal and sample topography information are recorded simultaneously. Measurements are always made while scanning in the same direction, i.e., positive y (y+) direction in order to avoid any appearance of effects of the backlash of the motors. While returning in y (y-), PID control is turned off and the z-piezo is fully retracted for safety. The next x-step is also made with the z-piezo fully retracted and PID control off. Then the PID control is turned on again and the next scan along y+ is performed. In the resulting SHPM images, the (x+) direction is from left to right and (y+) direction is from top to bottom. Note that it is the sample which is moved during the scan while the probe remains stationary.

In our setup, apart from this 'regulation mode', scanning at a fixed height from the sample surface in 'flyover mode', without such feedback, is also possible. Typically, in a particular measurement, a first scan is performed in the regulation mode, with the probe close to the sample surface. Next, the probe is retracted to the desired height from the sample surface and a scan is performed without topographic feedback for regulation of the extension of the z-piezo stack. Instead, the extension is determined from the topography information recorded earlier in the regulation mode. In this way, a set of scans can be performed at various heights from the sample surface. Note that this flyover mode of scanning allows to follow the topography and tilt of the sample. This is unlike standard lift-off mode scanning where the probe is lifted to a fixed height above the sample surface and scanning is performed in a flat plane without taking into account the sample topography and tilt. Standard lift-off mode scanning in a flat horizontal or tilted plane is also possible in our microscope. It is possible to perform scans at a height of up to 35 mm from the sample surface (range of the Z-motor).

B(z) profiles

Apart from performing scanning measurements, the system is also capable of, and has been extensively used for, fast recording of local B(z) profiles. For these measurements, the probe is approached close to the sample surface at the desired (x, y) location. Then the sample is retracted to a set probe-sample distance while the Hall signal is recorded at specified z intervals. In this way, local B(z) profiles at different desired (x, y) locations can be obtained.

III. MEASUREMENTS

The SHPM setup has been extensively used for quantitative characterization of a variety of magnetic structures, both hard and soft magnetic materials, topographically as well as thermomagnetically patterned micro-magnetic structures, and bulk permanent magnets.

Discussed in the following sections are some of the experimental results highlighting the various capabilities of the microscope.

A. 1D array of topographically patterned soft micro-magnets structures positioned above a bulk permanent magnet. Details of the preparation and characterization of the samples can be found in Ref. [START_REF] Roy | Fabrication and characterization of polymer membranes with integrated arrays of high performance micro-magnets[END_REF]. The scan resolution in Figures 5(a of the z-component of the stray magnetic field produced above such structures made with isotropic hard magnetic spheres [START_REF] Dempsey | Micro-magnetic imprinting of high field gradient magnetic flux sources[END_REF].

C. Bulk magnets

A SHPM system is a potentially attractive tool for characterization of bulk magnets, as on one hand it can detect large magnetic fields, and on the other hand, it can reveal local inhomogeneities that could affect the magnet performance. Our setup has been utilized

for characterization of a number of mm sized bulk magnets. Figure 8 shows a summary of measurements on a cubic magnet (5 mm side length). SHPM images in Figures 8(a)

and 8(b) show the B z distribution across a 2 mm × 2 mm area centered near the center of the magnet, at scan height h = 100 µm and 500 µm (± 1 µm), respectively. Note that these mentioned heights are the heights of the probe from the sample surface at the location (x, y) = (0, 0), i.e., the top-left of the image. The actual height at each point would vary according to the roughness of the magnet surface as in these sets the scans are performed in smooth horizontal planes (in flyover mode, using an artificial flat topography for reference).

This ensures surface corrugations do not affect the measured B z profiles. The scan resolution is 2.5 µm (image area 800 × 800 pixel). Figure 8(a) demonstrates inhomogeneities in the magnetic field distribution across the surface of the bulk magnet. The B z value is seen to vary in the large range of 360 mT < B z < 540 mT over the image area. Moreover, pockets of low B z (e.g., cf. dark blue region labelled 2 in Figure 8(a)) are observed. These observations are useful because these bulk magnets are used for a variety of purposes, e.g., to magnetize the soft magnetic microstructures discussed in Section III A. In view of the inhomogeneous B z distribution over the magnet surface, when such bulk magnets are used to characterize flux concentration in soft micro-magnets, it becomes important to take care in choosing the region of the surface onto which the micro-magnet array is placed. However, the image in Figure 8(b) at a larger height of 500 µm from the magnet surface shows that the B z distribution is largely smoothened out at this height from the magnet surface. This means at these heights, the issue of magnetic field inhomogeneity is not as crucial as in the earlier case.

This is further clarified from local B z (z) profiles at different locations on the magnet. With increasing z, the area over which the averaging occurs increases, and hence local 'anomalies' are averaged out. In general at a distance of ∼ 1-2 mm the B z distribution becomes homogeneous over an image size corresponding to the height, though there is still noticeable difference between the local B z values even at heights as large as 5 mm (cf. Figure 8). Beyond z ∼ 10 mm, all the curves merge into a single one, indicating uniform B z (z) across the magnet surface over the 2 mm × 2 mm area at these heights. Note that in both images in Figures 8(a) and 8(b), there is an overall gradient in the B z distribution along the y-direction (top-to-bottom of the images). This is attributed to the fact that the Hall probe is kept at a tilt of α ∼ 10 • w.r.t. the sample surface, the angle being in the yz-plane. Hence the effective measured component of B z changes as the probe is scanned along the y-axis.

Another factor contributing to this gradient could be that the scan area of 2 × 2 mm was not exactly centered at the middle of the magnet. Note that there is no such gradient along the x-axis. This issue is particularly noticeable in this case due to the intense fields and large area of the magnet. In usual scenarios of smaller sized micro-magnet samples this can be ignored, as seen in Sections III A and III B.

IV. CONCLUSION

We have developed a Scanning Hall Probe Microscope capable of performing high resolution scans (step size of 0.1 µm) over large areas (few mm, limited by sample topography and tilt). New generation of mesa etched and precision diced Hall sensors enable approaching close to the sample surface (< 2 µm) resulting in images with highly resolved features (spatial resolution of magnetic field distribution limited by the size of Hall probes, 1 µm).

These Hall probes allow detection of large magnetic fields (∼ 1 T) with high field resolution (100 µT, limited by system noise during scanning). The AFM-like tuning fork based force detection technique enables imaging in regulation mode (probe in contact with the sample surface) and the microscope can be also used for imaging in flyover mode at variable heights from the sample surface (up to 35 mm). Apart from scanning, the system is capable of recording local B(z) profiles. A single instrument (the Zurich Instruments HF2LI Lockin amplifier) serves as the primary electronics component of the microscope which significantly improves system performance. New programs have been developed for reliable operation and several enhanced functionalities such as precise sensor calibration, more efficient scanning in sets, fast acquisition of B(z) profiles etc. The simple design, three large range, precise stepper motors and a simple piezo stack combined with a reliable tuning fork based height control makes this microscope a very versatile tool for the control and inspection of devices and materials. The setup is being extensively used to measure stray field patterns of micro-magnet arrays as well as for characterization of bulk magnets.

FIG. 1 .

 1 FIG. 1. (a) Schematic of the microscope. (b) Optical image of the microscope showing the probe and sample stages. (c) The sample stage showing the Micos motors.

  5 mm (cf. results of our measurements on bulk magnets in Section III C showing 2 mm × 2 mm images). For comparison, minimum step size and maximum possible scan range at room temperature of some other reported stepper motor based long-range SHPM systems are as follows: 1.25 µm and 25 mm × 25 mm in Ref. 8, 0.2 µm and 1 cm × 4 cm in Ref. 9, 5 µm and 7 mm × 25 mm in Ref. 10, and 0.1 µm and 1 mm × 1 mm in Ref.11. All three motors in our setup are characterized by very reliable repeatability (negligible backlash error). The linear stages are controlled by a "SMC Corvus PCI"controller [22] installed in the PC. The system has a fast response, allowing stepping at 50 ms intervals. At this stepping interval, a typical scan (200 × 200 pixel) takes about two hours with additional overheads (cf. Section II C 2 for details on the scan protocol). In comparison, in the SHPM system described in Ref. 9, a 200 × 200 pixel scan takes several hours, and in Ref. 8 a 60 µm × 60 µm scan with pixel size 1.25 µm takes about 90 minutes.

FIG. 3 .

 3 FIG. 3. (a) and (b) Frequency response of a tuning fork glued to a probe plate, without and with a probe mounted, respectively. (c)Approach curve of a tuning fork with a probe mounted on it (cf. text for details).

FIG. 4 .

 4 FIG. 4. (a) Optical image of a Hall bar. (b) Optical image of a Hall bar showing the region of the Hall bar marked by the red rectangle in (a). The Hall cross is indicated by the red arrow. (c) SEM image of the Hall bar showing the Hall cross and leads. The Hall cross is indicated by the red lines. The distances of the Hall cross from the edges of the Hall bar are labelled as L 1 and L 2 . (d) Typical calibration curve of a Hall probe showing B vs. R Hall data and linear fit of the data.(e) and (f) Side and front views, respectively, of the probe, showing the angles α and β the probe makes with the sample surface.

  (a) shows an optical image of a typical Hall bar. The optical image in Figure 4(b) shows a close-up of the Hall bar with the Hall cross indicated by the red arrow.Figure 4(c) shows an SEM image of the Hall bar showing the Hall cross and contact leads. The Hall cross is identified by the red lines drawn on the image. The four edges of the Hall cross are marked as 1, 2, 3 and 4. The distances of the Hall cross from the edges of the Hall bar are shown as L 1 and L 2 . Hall probes of size 1-5 µm have been prepared. The particular probe shown here is of size ∼ 4 µm. Resistance values through the probes were found to be ∼ 15-17 kΩ for the smaller probes and 5-6 kΩ for the larger probes.

Figure 4 (

 4 Figure 4(d) shows typical calibration data for one of the Hall probes. The y-axis shows the

For proper interpretation of

  experimental data, it is important to reliably estimate the distance between the Hall cross and the sample surface. To avoid contact of the wires micro-bonded to the Hall probe with the sample surface and for safety of the probe against damage by hitting the sample surface, the probe is kept at a slight inclination w.r.t. the sample stage. Figures4(e) and 4(f) show side and front views, respectively, of a probe mounted on a tuning fork with a sample below, indicating the angles α and β the probe makes with the sample surface. Note that the Hall cross is not exactly at the edge of the bar but at a distance from it. In Figure4(c) the distances of the Hall cross from the Hall bar edges are indicated as L 1 and L 2 . For the smallest probes (1 µm size), L 1 = L 2 = 8 µm. For the largest probes (5 µm size), L 1 = L 2 = 15 µm. The angle α ensures that the probe does not hit the sample surface and the angle β ensures that the edge of the Hall bar closer to the Hall cross touches the sample surface. The angle α is generally set at ∼ 9-10 •

FIG. 5 .

 5 FIG. 5. (a)-(g) SHPM showing the B z distribution in a 1D array of topographically patterned Fe-Co soft micro-magnet pillar structures embedded in a polymer (cf. text for details). The image in (a) is obtained in the regulation mode and those in (b)-(g) are obtained in the flyover mode (cf. text for details). The height of the probe from the sample surface corresponding to each image is indicated in the respective panels. The inset in (g) shows the common scale bar. (h) Image showing topography of the sample. It is obtained simultaneously along with the SHPM image during the regulation mode scan.

Figures 5 (

 5 Figures 5(a)-(g) show a set of SHPM images showing the z-component of the stray field pattern (B z ) in a 1D array of topographically patterned Fe-Co soft micro-magnet pillar

FIG. 6 .

 6 FIG. 6. (a) and (b) B z (y) profiles at different scan heights along the red lines shown in Figure 5(a), one profile along a region with high field gradients and one along a region with low field gradients. (c) B z (z) profiles at different locations on the image as indicated by numbers in 5(a).

  )-(g) is 0.5 µm and the scan area = 400 × 200 pixel = 200 µm × 100 µm. The B z color scale in the images is indicated by the common scale bar in the inset of Figure 5(g). Each micro-pillar is 10 µm wide and the separation is 30 µm (center-to-center). The pillars are 40 µm tall. For these measurements, the sample is placed on top of a bulk magnet which produces a field of ∼ 320 mT normal to the sample plane at the sample surface. The first image in Figure 5(a) was obtained in the regulation mode with the probe-sample distance estimated at h ∼ 2 µm (± 0.5 µm). Then successive images were obtained at different probe-sample distances in the flyover mode up to h ∼ 68 µm. The h corresponding to each image is indicated in the respective panels. The topography recorded in the regulation mode is shown in Figure 5(h). The Fe-Co pillars are coated with a layer of PDMS that smoothens the large topographical variations in the sample due to the micro-pillars and hence enables safe scanning of the Hall probe over the sample surface. In Figure 5(h) modulations in the PDMS layer due to the underlying micro-pillars, along with an overall tilt in the sample plane w.r.t. the scanning plane, can be observed. From Figures 5(a)-(g), strong field modulation, with large flux concentration at the micro-pillars can be clearly observed, with B z as high as 713 mT measured 2 µm above the sample surface at the position of the pillars (cf. Figure 5(a) and scale bar in inset of Figure 5(g)). In Figure 5(a), a corresponding decrease in flux density in the vicinity of the pillars can also be observed. With increasing probe-sample distance h, the field modulation due to the micro-pillars can be seen to reduce. Finally, at h ∼ 68 µm, modulations completely disappear, and a uniform field distribution across the image, corresponding to the underlying bulk magnet, is observed. These observations are further elaborated in Figure 6. Figures 6(a) and 6(b) show B z (y) profiles obtained from the set of images shown in Figure 5 along the two red lines drawn over the image in Figure 5(a), through the center of a pillar and in between two pillars, respectively.The height h corresponding to each B z (y) profile is indicated in both panels. The strong flux concentration at the micro-pillars is very clear from Figure 6(a). Figure 6(b) clearly shows the dip in B z close to the pillars. Further, B z (z) profiles were recorded at various locations of the image area as indicated by numbers in Figure 5(a). The B z (z) profiles are shown in Figure 6(c). From Figure 6(c), the strong flux concentration at the micro-pillars (> 700 mT near the sample surface at positions 2 and 3, on top of micropillars) over and above the field due to the bulk magnet (∼ 320 mT near the sample surface at positions 1 and 7, far from the micro-pillars) is very clear. Further, the B z (z) profiles show that the field modulation effect of the pillars completely disappears beyond z ∼ 0.1 mm, where all the profiles merge into the profile expected for the bulk magnet only. B. Hard magnetic powder based micro-flux sources We have characterized hard magnetic powder (NdFeB) based micro-flux sources prepared using the micro-magnetic imprinting (µMI) technique [35]. More details on the µMI tech-nique and preparation and characterization of micro-flux sources by this method can be found in Ref. 35. The images in Figures 7(a)-(c) show the z-component of the stray field distribution (B z ) produced at heights h = (a) 5 µm, (b) 20 µm, and (c) 50 µm (± 1 µm) above a µMI structure of an array of 100 µm × 100 µm squares of NdFeB powder (average NdFeB particle size = 5 µm). For each image, the scan resolution = 2.5 µm and the scan area = 160 × 400 pixel = 400 µm × 1000 µm.

Figure 7 (

 7 d) shows B z (x) profiles across the dashed lines shown on each image. Close to the sample surface, the SHPM technique is able to clearly resolve the well-defined individual square structures and large stray field values are detected (cf. image in Figure 7(a) at h = 5 µm and peak-to-peak value of ∼ 40 mT in the corresponding B z (x) profile in Figure 7(d)). Further, smaller agglomerations of powders across the sample surface are also clearly observed in the image. With increasing scan height, as expected, the individual structures become less resolved and the detected stray fields also decrease. The measurements have been shown to agree very well with simulations

FIG. 7 .

 7 FIG. 7. (a)-(c) SHPM images showing the B z distribution above a µMI structure of an array of 100 µm × 100 µm squares of NdFeB powder at height h = 5, 10 and 20 µm, respectively. The individual color bars are shown alongside each image. (d) B z (x) profiles across the dashed lines shown on each image.

FIG. 8 .

 8 FIG. 8. (a) and (b) SHPM images showing the B z distribution above a 5 mm cubic bulk magnet at scan height 100 and 500 µm, respectively. The individual color bars are shown alongside each image. (c) B z (z) profiles at different locations on the image as indicated by numbers in (a).

Figure 8 (

 8 Figure 8(c) shows B z (z) profiles obtained at five different locations on the image as indicated by the numbers in Figure 8(a). Interestingly, while the profiles at the corners (locations 1, 3 and 4) and center of the image (location 5) show expected monotonic decrease in B z with increasing z, in the profile at the center of the low B z region (location 2), the B z value initially increases with increasing z. This continues till z ∼ 0.4 mm, beyond which the curve reverses and follows the same trend as the others. With increasing z, the difference in local B z values at different locations decreases. This is because the detected 'local' field really is an average of field lines emanating from the area surrounding a given location.

  johnson ∼ 1.424 nV/ √ Hz. For excitation current of 10 µA, this gives a lowest possible noise in B z of ∆B z ∼ 1 µT/ √ Hz, or a magnetic field sensitivity of 10 µT (at measurement frequency 522 Hz and time constant 10 ms). However, the magnetic field resolution in our measurements is ∼ 100 µT. Thus additional noise sources other than the Johnson noise are at the origin of this field resolution. A noise floor of ∼ 16 nV/ √ Hz is measured when the Hall probe is not excited. It increases to ∼ 18 nV/ √ Hz when the probe is excited. The measured noise level is thus essentially due to the wiring of the probe, capturing ambient electrical noise, as the microscope operates in an unshielded environment. The field sensitivity in our system is in between values reported for SHPM systems operating at room temperature (300 K): 0.08

	µT/ √	Hz in Ref. 8 (InSb Hall probes), 350 µT/ √	Hz in Ref. 11 (GaAs/AlGaAs Hall probes).

C. Measurement protocol

1. Estimate of distance of probe from sample surface

ACKNOWLEDGMENTS This work was partially funded by the ANR "MIME"project (ANR-11-BSV5-0101). The work of R. B. G. K. is partially supported by the "Chaire Initiative Universitaire Alpes (IUA)"of the University Grenoble Alpes. The authors would like to thank T. Devillers, D Le Roy, F. Dumas-Bouchiat, H. Marelli-Mathevon, A. Dias, and D. Givord for providing samples and for useful discussions. G. S. would like to thank the electronic workshop at Institut Néel for technical support and S. Le-Denmat for his inputs. The authors acknowledge the contributions of D. J. Hykel, M. Kustov and P. Laczkowski in development of the earlier generation of the microscope. The 2DEG wafers used for probe fabrication were provided by Picogiga (Soitec). Probe patterning was carried out at the Nanofab cleanroom facility of Institut Néel.