Quasiparticle Diffusion and Vortex Detection Models for SNSPD

INOM

EXAMENSARBETE TEKNISK FYSIK, AVANCERAD NIVÅ, 30 HP

STOCKHOLM SVERIGE 2020,

Quasiparticle Diffusion and Vortex Detection Models for SNSPD

YANI SHEN

KTH

SKOLAN FÖR TEKNIKVETENSKAP

Abstract

This thesis presents a theoretical study of the photon detection mechanism in su- perconducting nanowire single-photon detectors, with focus on modelling the process of normal-conducting domain generation. Investigation of the subject led to some mod- ifications on the existing model used in recent articles. Simulations showed that the relaxation of supercurrent back to a homogeneous state after photon absorption happens faster for photon absorption at the center of the superconducting strip than absorption near the strip edge. It was also found that in case of photon absorption at the strip center the single vortex entry barrier has two local maxima, one of which disappeared with increasing time. For photon absorption near the strip edge the entry barrier instead has only one maximum.

Den här uppsatsen presenterar en teoretisk studie av fotondetektionsmekanismen i supraledande nanotråd singel-foton detektorer, med fokus på modellering av processen för generation av ett normalledande område. Undersökning av ämnet ledde till några modifikationer av nuvarande modell som används i nyligen publicerade artiklar. Simu- leringar visade att relaxationen av superström tillbaks till ett homogent tillstånd efter fotonabsorption händer snabbare för fotonabsorption i mitten av den supraledande trå- den än för absorption nära trådkanten. Det visades också att i fallet med fotonabsorption i mitten av tråden har singel vortex inträngningsbarriären två lokala maxima, där det ena försvann med tiden. Inträngningsbarriären har bara ett maximum för fotonabsorption nära trådkanten.

Acknowledgements

This thesis is the result of six months of full-time studies, for the degree of Master of Science in Engineering Physics, at the Department of Physics at the Royal Institute of Technology, Sweden.

I would like to thank my supervisor Prof. Mats Wallin for giving me the opportunity to work on this project and for providing excellent guidance and support.

I would also like to thank my family and my friends for their support.

Contents

1 Introduction 1

2 Superconductivity 3

3 Model of normal-conducting domain generation 5

3.1 Quasiparticle diffusion . . . 5

3.2 Supercurrent redistribution . . . 7

3.2.1 Model of supercurrent redistribution . . . 7

3.2.2 Implementation of supercurrent redistribution . . . 8

3.3 Triggering of voltage pulse . . . 10

3.3.1 Single vortex crossing . . . 11

4 Simulations and discussions 14 4.1 Photon absorption at strip center . . . 14

4.1.1 TaN superconducting strip . . . 14

4.1.2 NbN superconducting strip . . . 23

4.2 Photon absorption near strip edge . . . 31

4.2.1 TaN superconducting strip . . . 31

4.2.2 NbN superconducting strip . . . 38

5 Summary 46

Bibliography 48

Chapter 1 Introduction

A single-photon sensitive optical detection device generates an electrical signal upon ab- sorption of a photon. There are lots of performance metrics that quantifies the use and limitations of such devices, among which the most important ones are detection efficiency, dark count rate, time jitter, and recovery time [8],[11]. The detection efficiency measures the probability of electrical signal registration if a photon is absorbed by the detector.

The dark count rate measures the probability of electrical signal registration without photon absorption, which could happen due to for example electrical noise, thermal fluc- tuations, or a too high bias current [2],[8]. The time jitter, or time uncertainty measures the variation in time interval between photon absorption and signal generation. The recovery time, or dead time, is the time it takes for the resistive state to relax back into the superconducting state to prepare for next photon registration. A shorter recovery time implies a higher photon count rate. An effective optical detection device should have high detection rate, low dark count rate, small time jitter, and short recovery time.

There is a variety of single-photon detection devices [8]. The most promising ones are semi-conducting avalanche photodiodes (APDs) and superconducting nanowire single- photon detectors (SNSPDs) [15]. One drawback of SNSPD is its relatively low operating temperature of few Kelvins, whereas most APDs have a much higher operating tem- perature of few hundred Kelvins, which makes APD relatively inexpensive compared with SNSPD [15]. Another drawback of general SNSPDs is their insensitivity of photon count number, i.e., they are able to register either zero photon or "one or more than one photons". However, there is a special SNSPD configuration, the so-called photon- number-resolving SNSPD (PNR-SNSPD), that deals with this problem [11]. Overall, SNSPD compares well with APD, due its high performance in all aforementioned perfor- mance metric aspects [11],[15]. Reference [15] reported a NbTiN SNSPD with a photon detection rate over 150 MHz and a dark count rate below 130 Hz. The time jitter could be reduced to 14.80 ps while keeping the overall detection efficiency higher than 75%.

SNSPD has often been the detector of choice in many demanding applications such as quantum key distribution, which is one type of quantum cryptography, where a low dark count rate and a small time jitter is required [13]. Another example is application in space-to-ground communication: an optical receiver based on NbN SNSPD arrays was under development for the NASA Lunar Laser Communications Demonstration program [7].

SNSPD consists of a thin, narrow, and long superconducting nanowire, which is main- tained at a temperature T below its superconducting critical temperature T_c, and biased with a constant direct current Ibias below its superconducting critical current Ic. Assume a photon with energy E is absorbed by the current-biased detector. If E and I_bias are in correct combination, a voltage pulse is triggered which is then registered as a photon count. For a given bias current, there is a minimum photon energy, in other words, a max- imum photon wavelength, above which the detection mechanism can not be triggered.

Also, the photon energy should not exceed the work function of the superconducting film.

SNSPD is typically sensitive at visible and infrared wavelengths [8],[11].

The photon detection process in SNSPD may be divided into three steps: photon ab- sorption, normal-conducting domain generation, and electrical signal output [4]. The first step deals with understanding and improvement of the probability of photons actually being absorbed by the detection device. It is usual to form the nanowire in a meander pattern to maximize the area of photon incidence. However, not all photons incident on the device will be absorbed, some will be reflected from or transmitted through the detector. One solution is to integrate the SNSPD device with an optical cavity and an anti-reflection coating [12]. The last step deals with the property of readout electron- ics. It has been shown that the electrical readout circuit can set a limit on the photon count rate and how this problem can be solved by improving the design of the electrical circuit [10]. The focus of this project is on the second step: normal-conducting domain generation. This step can in turn be divided into three parts: quasiparticle diffusion, su- percurrent redistribution, and triggering of voltage pulse [5]. All steps will be explained further in chapter 3. The model of this process is formulated by following recent articles of the relevant subject.

Despite intense experimental and theoretical studies detailed understanding of the photon detection mechanism of SNSPD, especially the mechanism responsible for the triggering of normal-conducting phase, remains incomplete. The goal is to find a the- ory that gives detailed understanding of the photon detection mechanism, and thereby help to further improve the performance of SNSPD by optimizing the detector material, detector geometry etc. The aim of this project is to gain more understanding of the normal domain generation process and thereby take a step towards this goal. During the investigations, a more suitable approach for calculating supercurrent redistribution was found. The contribution of vortex core energy, which has been neglected in some earlier published articles [2],[3],[4],[5],[9], is included here.

This thesis is structured in the following way. Chapter 2 introduces the SNSPD de- tector geometry considered in the model and summarizes results for superconductivity that are relevant for the detector model. Chapter 3 describes the model and discusses the normal-conducting domain generation. Chapter 4 presents the results from numerical calculations. Chapter 5 gives a summary and conclusions.

Chapter 2

Superconductivity

Here we state relevant formulas in superconductivity that enters the numerical simula- tions of the normal domain generation model. All formulas are given in SI units. For derivations see Ref. [14].

Figure 2.1 below shows the geometry and coordinate system used for a straight wire segment of SNSPD. Following the typical length scale of such devices [5],[6], we set here the strip thickness d to 5 nm, width w = 100 nm, and length L = 1000 nm. The wire is sufficiently thin (d  w, L) so we could restrict the calculations to two dimensions. We set the x-axis to point along the direction of the applied bias current Ibias, and the y-axis transverse to the direction of I_bias.

The superconducting order parameter Ψ = |Ψ|e^iθ is the wave function averaged over large atomic distances so that variations on atomic distances is ignored. In a uniform state we have |Ψ|² = n_se,0 = N_se,0/V , where n_se,0 is the superconducting electron den- sity at zero bias current, Nse,0 the number of Cooper pairs, and V the sample volume.

The two dimensional density n^2D_se,0 is defined as N_se,0 = n_se,0V = n^2D_se,0×area, and with V = d×area this gives

n^2D_se,0 = n_se,0d, (2.1)

with n_se,0 given by

n_se,0 = me

µ₀q²λ², (2.2)

where m_eand q are the electron mass and charge, respectively. λ is the London penetra- tion depth, which is a temperature-dependent parameter given by

λ(T ) ≈ λ(T = 0)1 − τ⁴−1/2

, (2.3)

where τ = T /T_c is the reduced temperature.

The pearl length, or the effective magnetic penetration depth, is given by

Λ = 2λ²/d. (2.4)

For dirty superconductors, the electric mean free path is much less than the coherence length, resulting in a much larger value of magnetic penetration depth [1]. In the dirty limit, the London penetration depth λ is related to the Ginzburg-Landau penetration depth λ_GL by

λ(T ) = λGL

 ∆(τ )

∆ tanh ∆(τ ) 2k_BT

−1/2

, (2.5)

where ∆ is the superconducting gap given by

∆(T ) = ∆(0)1 − τ²0.5

1 + τ²0.3

. (2.6)

The zero-temperature superconducting gap is ∆(0) = αk_BT_c, where the BCS value α = 1.764 is assumed for TaN, and α = 2 for NbN [5].

The coherence length ξ is also temperature-dependent [5]

ξ(τ ) = ξ(0)(1 − τ )^−0.5(1 + τ )^−0.25. (2.7)

Figure 2.1: SNSPD geometry and the coordinate system used for a straight wire segment.

Typically d  w  L, so we could restrict the calculations to 2D. The bias current Ibias

is set to point along the x-axis.

Chapter 3

Model of normal-conducting domain generation

This chapter explains the three parts of the normal-conducting domain generation model, which are quasiparticle diffusion, supercurrent redistribution, and triggering of voltage pulse [5].

3.1 Quasiparticle diffusion

Following article [5], we start with quasiparticle diffusion as the first step after photon absorption. The idea is that the energy of a photon absorbed by the superconducting film surface will be transferred to one electron, which starts to diffuse away from the photon absorption site according to the classical diffusion equation

∂C_e(~r, t)

∂t = D_e∇²C_e(~r, t), (3.1)

where D_e is the diffusion coefficient of normal electrons and C_e(~r, t) is the probability density to find the excited electron at position ~r at time t.

The excited electron undergoes inelastic scatterings while it diffuses. The photon energy sensitive to SNSPD is usually several eV, whereas the typical superconducting bandgap is on the scale of meV. Therefore, a certain fraction of these scattering events will break up Cooper pairs and thereby cause the creation of quasiparticles. This fraction is represented by the conversion efficiency ζ, which has to be determined experimentally.

Here ζ is set to the value 0.25.

The quasiparticles themselves also undergo diffusion within the superconducting strip according to the classical diffusion equation, but with two extra terms

∂C_qp(~r, t)

∂t = D_qp∇²C_qp(~r, t) − C_qp(~r, t)

τ_r + ζE

∆τ_qp

n^2D_se,0− C_qp(~r, t) n^2D_se,0

 exp



− t τ_qp



C_e(~r, t).

(3.2)

The first term on the right hand side of Eq. (3.2) comes from the usual diffusion equation with quasiparticle diffusion coefficient D_qp. The second term gives a negative contribution to the quasiparticle probability density Cqp(~r, t) due to recombination of quasiparticles to Cooper pairs on a time scale τ_r. The last term accounts for the aforementioned cre- ation of quasiparticles due to scatterings of the excited electron. The excitation energy of the electron is assumed to decay exponentially with a time scale τqp. E is the ab- sorbed photon energy and ∆ is the superconducting energy gap. The multiplicative term n^2D_se,0− C_qp(~r, t)
/n^2D_se,0 accounts for the saturation of the number of quasiparticles. This term is not included in Ref. [5] but is taken from Ref. [9]. All material parameters used are listed in Tab. 3.1, taken directly from Ref. [5]. Note that since two quasiparticles are required to form one Cooper pair, the recombination term in Eq. (3.2) should be proportional to Cqp(~r, t)² instead of Cqp(~r, t). But as motivated in Ref. [5], because the recombination time scale τ_r is much larger than the energy decay time scale τ_qp that appears in the same equation, the recombination term will only have a minor effect so the model could be kept linear.

Equations (3.1) and (3.2) form a set of coupled differential equations. In order to solve these, one also need initial conditions and boundary conditions. A possible initial condition is to set C_e(~r, t₀) in the form of a narrow Gaussian function

C_e(~r, t₀) = 1

√4πD_et₀ exp



−(x − xa)²+ (y − ya)² 4D_et₀



, (3.3)

where t₀ > 0 is the initial time, x_a and y_a are the coordinates of the photon absorption site. This condition ensures initial electron normalization for an infinite 2D strip.

Another possible choice is to set C_e(~r, t₀) = 1/∆A on the photon absorption site, and zero on all other sites. Here ∆A = A/N , where A is the total area of the superconducting strip and N is the number of lattice sites when discretizing the strip into 2D grids. This ensures the initial normalization R C_e(~r, t₀) dxdy → P

iC_e,i∆A = 1, where i is the index over lattice sites. Initially there are no quasiparticles so C_qp(~r, t₀) = 0 in the whole strip.

The following boundary conditions are used

C(x = |L/2|, y, t) = 0, (3.4)

∂C(x, y, t)

∂y    y=|w/2|

= 0. (3.5)

The zero Dirichlet boundary condition represents particles flowing out of the segment, and the zero Neumann boundary condition means that particles can not pass through the strip side-walls. Equations (3.4) and (3.5) apply for both Ce and Cqp.

After solving the coupled diffusion equations, the result C_qp(~r, t) can be used in the next step of the model, which is the calculation of supercurrent redistribution.

∆ D_qp D_e ξ λ_GL τ_r τ_qp N₀ ζ meV nm²/ps nm²/ps nm nm ps ps nm⁻³eV⁻¹

TaN 1.3 8.2 60 5.3 520 1000 1.6 48 0.25

NbN 2.3 7.1 52 4.3 430 1000 1.6 51 0.25

Table 3.1: Material parameters used in the simulation for T = 0.05T_c. The density of states is here defined such that E_cond = ¹₂V N (0)∆² is the total condensation energy of the superconducting electrons in a volume V .

3.2 Supercurrent redistribution

In the absence of quasiparticles the supercurrent in the straight wire segment will be uniformly distributed. After photon absorption, currents will be redistributed due to the inhomogeneous distribution of quasiparticles throughout the superconducting strip. In regions with high quasiparticle concentrations, superconductivity is suppressed causing an increase of the current outside these region in order to maintain current conservation.

In this project, the model for calculating the current redistribution is different from the one used in articles [4],[5],[9]. The main difference is that we consider not only current continuity, but also take into account the minimization of current kinetic energy. The reason is that there are infinitely many current configurations that obey current continu- ity, which means that the condition of current conservation is necessary but not sufficient.

3.2.1 Model of supercurrent redistribution

The superconducting order parameter can be written as Ψ(~r) = |Ψ(~r)|e^iθ(~r), with the density of superconducting electrons being given by n^2D_se (~r) = |Ψ(~r)|². The supercurrent is given by

J =~ 1

2m(Ψ^∗P Ψ + Ψ( ~~ P Ψ)^∗) = n^2D_se

m (~∇θ − q ~A), (3.6)

where m = 2meis the mass of one Cooper pair and ~P = −i~∇−q ~A is the gauge invariant kinetic momentum operator. The total kinetic energy is given by

H = Z

Ψ^∗ P~² 2mΨ =

Z 1

2m|(−i~∇ − q ~A)Ψ|² = ~² 2m

Z

|∇|Ψ||²+

Z 1

2mn^2D_se |(~∇θ − q ~A)|², (3.7) where the first equality is by definition. To get the second equality use the fact that ~P is a Hermitian operator. The quasiparticle density C_qp(~r) at any time t determines the superconducting electron density n^2D_se (~r) at time t, see Eq. (3.14). Therefore, ∇|Ψ| is a fixed function regarding to supercurrent redistribution. In other words, the integral in Eq. (3.7) involving this term is unrelated to supercurrent redistribution so we can simply drop it. Rewriting the integrand of the second integral in Eq. (3.7) using Eq. (3.6), we obtain

H =

Z m

2n^2D_se | ~J |². (3.8)

So the kinetic energy of the supercurrent can be written as H =

Z K(~r)

2 | ~J (~r)|², (3.9)

where K(~r) = m/n^2D_se (~r). The discretized version of Eq. (3.9) on a square grid is

H =X

l

Kl

2 (J_x,l² + J_y,l² ), (3.10)

where l is the index over lattice sites.

A bias current I_bias = wJ_bias is imposed over the ends of the straight wire and is assumed to stay uniform far away from the photon absorption site. In the interior of the wire current redistribution is to be calculated by imposing the condition of local current conservation ∇ · ~J = 0 and by minimizing the kinetic energy H in Eq. (3.10).

3.2.2 Implementation of supercurrent redistribution

At the ends of the wire far enough from the photon absorption site, the current density is assumed to be uniform and is taken to be an input parameter. From this input the current distribution in the whole wire has to be calculated for a given quasiparticle dis- tribution. Distinguish from the wire ends at x = ±L/2 and the interior with |x| < L/2.

Assume as a boundary condition on the current a uniform distribution at the ends, given by Jx = J₀, J_y = 0, where J₀ is the input parameter corresponding to the strength of the bias current density. At all interior links the current distribution is calculated by the following simple steps.

1. As initial condition, set J_x = J₀, J_y = 0 also on the interior lattice links. This state obeys current conservation but does not minimize the kinetic energy H if the density of superconducting electrons n^2D_se is not uniform.

2. Relax the current distribution until the minimum of H is reached by the following simple iteration. Look at all interior elementary plaquettes in the lattice. See Fig. (3.1).

The idea is to add a directed plaquette loop current ∆ to minimize H locally. Using Eq. (3.10), the kinetic energy can be written as

H = 1

2K_x(i, j)(J_x(i, j) − ∆)²+ K_y(i + 1, j)(J_y(i + 1, j) − ∆)²

+K_x(i, j + 1)(J_x(i, j + 1) + ∆)²+ K_y(i, j)(J_y(i, j) + ∆)² . (3.11)

Figure 3.1: A directed plaquette loop current ∆ used to minimize the supercurrent kinetic energy H. Note that circulation of ∆ obeys current continuity.

The value of ∆ is found by setting ^dH_d∆ = 0 which gives

∆ = [Kx(i, j)Jx(i, j) + Ky(i + 1, j)Jy(i + 1, j)

−K_x(i, j + 1)J_x(i, j + 1) − K_y(i, j)J_y(i, j)] /K, (3.12)

where K = K_x(i, j) + K_y(i + 1, j) + K_x(i, j + 1) + K_y(i, j). Since the only changes are circulating loop currents obeying continuity, i.e., as much in and out at each lattice site, the current remains continuous throughout the calculation. The energy H will relax monotonically towards its lowest value that gives the desired current distribution. Since this is a local relaxation method it has to be iterated by sequentially sweeping through the lattice until the increments ∆ are small enough. After each iteration, the interior link current densities are updated according to Fig. (3.1). The end link values are held fixed throughout the calculation. Next the local coupling constants K_x, K_y have to be specified. They are here related to the average quasiparticle density on the lattice sites connected by the links:

K_x(i, j) = 2m

n^2D_se (i, j) + n^2D_se (i + 1, j), K_y(i, j) = 2m

n^2D_se (i, j) + n^2D_se (i, j + 1). (3.13)

As in Ref. [5], it is assumed that the increase of the local quasiparticle concentration equals the decrease of the local superconducting electron concentration, i.e.,

n^2D_se (~r, t) = n^2D_se,0− C_qp(~r, t), (3.14) where n^2D_se,0 is the 2D superconducting electron density at zero bias current, given by Eqs.

(2.1) and (2.2). This is a linear model. A more accurate model was considered in Ref.

[4], where the fact that high drift velocities have a pair-breaking effect that reduces n^2D_se was included

n^2D_se → n^2D_se

 1 − 1

3

 J Jc

2

. (3.15)

Here the current has to be solved self-consistently by iteration.

The spatial and temporal evolution of the quasiparticle density Cqp(~r, t) obtained by solving the coupled differential equations can now be put into Eq. (3.14) to obtain the evolution of superconducting electron density, which in turn can be used to determine the supercurrent redistribution by following the numerical iteration process descried above.

Once we have n^2D_se(~r, t) and ~J (~r, t) at hand, we can move on to the last step of the normal domain generation model, which is the triggering of a voltage pulse.

3.3 Triggering of voltage pulse

Three different models of the voltage pulse trigger mechanism were considered in Ref. [5].

1. Hard-core model: assumes that the quasiparticle concentration is high enough to completely suppress superconducting electron density in a region around the photon absorption site. The current density in this region is assumed to be zero and we have a so called normal-conducting hot spot. The current inside the spot is diverted along the periphery of the hot spot into the still superconducting region outside as required by current continuity. If the current outside the spot exceeds the critical current the whole cross-section becomes normal-conducting leading to a voltage signal which is registered as a photon detection. For a given bias current, there is a minimum photon energy below which the detection can not be registered and vice versa.

2. QP model: this model is an improved version of the hard-core model, which did not take into account quasiparticles outside the normal-conducting hot spot. In this model, the normal region forms after when the critical current has dropped down to the applied current.

3. Single vortex crossing model: assumes that the photon absorption triggers a single magnetic vortex, which appears initially at the edge of the superconducting strip, to move across the width of the strip without the whole cross-section becoming normal- conducting. The vortex crossing leads to a 2π phase slip of the superconducting order parameter θ, which results in a burst of voltage that is detected and registered as a pho- ton count. The voltage is given by the Josephson phase relation

V = ~ 2e

dθ

dt. (3.16)

In Ref. [5], it was found that for a given photon energy, the vortex model leads to detection events at lower bias currents and the detection also occurs earlier in time com- pared to the other two models. Their result indicated a strong support of the entry and crossing of a single magnetic vortex as the primary photon register mechanism.

Another voltage pulse trigger model is creation and crossing of a vortex-antivortex pair instead of a single vortex. However, in Ref. [2], it was shown that the energy barrier for vortex-antivortex pair is higher than that for a single vortex. Therefore, we choose to focus on the single vortex crossing model. More detail of this model is described in the

next section.

3.3.1 Single vortex crossing

Here we present the existing formula of potential energy barrier for single vortex crossing taken from Ref. [5], but use a different normalization and add to it the contribution of vortex core energy.

In the case of a homogeneous current density, i.e., when J_bias = I_bias/w is constant across the strip length, then energy barrier prohibiting the single vortex entry and cross- ing is given by

U (y, I_bias, T )

0

= ln

 2w

πξ(T )cosπy w



−I_bias Ic,v

2(y + w/2)

exp(1)ξ(T ), (3.17) where I_c,v is the bias current for which the potential energy barrier for single vortex crossing vanishes in a straight wire. The potential energy on the LHS is scaled with the characteristic vortex energy ₀ given by

0 = Φ²₀

2πµ₀Λ = Φ²₀e²

4πm_en^2D_se,0, (3.18)

where the second equality comes from using formulas listed in chapter 2.

The first term on the RHS of Eq. (3.17) represents the vortex self energy and the second term is the work done by the bias current on the vortex [9]. The first term represents the Bean-Livingstone surface entry barrier due to the boundary condition of no current across the boundary. The vortex-vortex interaction is a logarithmic 2D Coulomb interaction, U (r) ∼ −q₁q₂ln(r/r₀), where q = ±1 is the vorticity corresponding to Coulomb charge. The boundary condition is therefore implemented as an image vortex outside the sample with opposite vorticity that attracts the vortex to the boundary. This surface barrier must be overcome for the vortex to enter the interior of the wire. The Lorentz/Magnus force from the bias current counteracts this surface barrier. This is why a higher bias current results in a lower vortex entry barrier U .

In case with a non-vanishing spatial-dependent quasiparticle density, the supercurrent density will no longer be uniform. Therefore we differentiate the potential energy U and integrate it again taking into account the position dependent quantities n_seand J_x. Also, the energy scale ₀ in Eq. (3.18) becomes position-dependent

₀(x, y) = Φ²₀e² 4πme

n^2D_se (x, y). (3.19)

The potential barrier is now given by U (x, y)

₀ = −π w

Z y

−^w

2+ξ

n^2D_se (x, y⁰)

n^2D_se,0 tan πy⁰ w



dy⁰− 2w I_c,vexp(1)ξ

Z y

−^w

2

n^2D_se (x, y⁰)

n^2D_se,0 j_x(x, y⁰)dy⁰. (3.20)

Note that the lower integration boundary of the first integral is at one coherence length away from the lower edge of the strip. In the sample edge regions the vortex is overlapping with its mirror image and the self energy diverges. This divergence is removed by setting the first term in Eq. (3.17) to zero within a distance ξ from the strip edges by the selection of integration limits. Physically the cutoff is due to the suppression of the supercurrent inside the vortex core which eliminates the divergence of the interaction at short distance.

In Ref. [5], the normalization condition U = 0 for y < ξ − w/2 and y > w/2 − ξ was used. This boundary condition is not correct and should be replaced by including the vortex core energy in U . The loss of condensation energy for making the vortex core normal increases the cost of inserting a vortex in the wire and including it can be expected to improve the existing model of the vortex energy barrier U . In the bulk of the film the total vortex core energy E_c is included, and at the sample edges the fraction of the vortex core inside the sample is included. If a vortex has penetrated a distance y < ξ from the edge, an energy Ecy/ξ should be included in U . This comes from a region of normal-conducting state of area A ≈ yξ at the sample edge with an energy corresponding to a fraction A/ξ² of the total core energy.

The condensation energy is the difference in energy between electrons in the normal- conducting state and the superconducting state. Therefore, the creation of a vortex with normal electrons in the vortex core increases the total energy with the condensation en- ergy of electrons within the vortex core, which we call the vortex core energy E_c. The condensation energy density is given by [14]

_cond = 1

2N (0)∆², (3.21)

where ∆ is the energy gap and N (0) is the density of states at the Fermi level listed in Tab. 3.1. Hence the vortex core energy is given by Ec = ¹₂V_coreN (0)∆², where Vcore is the core volume. Since the superconducting electron density n_se is non-uniform, we need to identify the position dependence of the condensation energy. Using that N (0)∆ ∼ n_se, the core energy for a vortex at position (x, y) with a cylindrical core volume Vcore= πξ²d can be estimated to be

E_c(x, y) ≈ 1

2πξ²dN (0)∆²n^2D_se (x, y)

n^2D_se,0 . (3.22)

Note that for vortex-antivortex pair, the core energy contribution is twice of E_c.

There is another important effect which a vortex core has on the potential energy:

suppression of superconducting electron density inside the vortex core leads to current crowding outside the core. The new current density should then be recalculated using the relaxation method described in Chapter 3.2. For practical purposes one can take n^2D_se to be a small number inside the vortex core. It is easy to implement this in the code for a vortex inserted by hand at some fixed position. However, to simulate vortex crossing motion through the strip taking this effect into account is computationally expensive.

This remains to be done as a future development of the calculation.

To summarize, the potential energy barrier for vortex crossing in case of non-uniform superconducting electron density used in simulations here is given by

U (x, y)

₀ = −π w

Z y

−^w

2+ξ

n^2D_se (x, y⁰)

n^2D_se,0 tan πy⁰ w



dy⁰− 2w I_c,vexp(1)ξ

Z y

−^w

2

n^2D_se (x, y⁰)

n^2D_se,0 j_x(x, y⁰)dy⁰+E_c

₀, (3.23) with Ecgiven by Eq. (3.22). The upper integration limit of the first integral is y = w/2−ξ and for the second integral the upper integration limit is w/2.

Chapter 4

Simulations and discussions

4.1 Photon absorption at strip center

4.1.1 TaN superconducting strip

The Matlab PDE-solver was used to solve the coupled diffusion equation system, Eqs.

(3.1) and (3.2), on a straight rectangular superconducting strip segment shown in Fig.

2.1. Simulation parameters are listed in Tab. 3.1. The time step was set to dt = 0.1 ps.

The incident photon was assumed to have a wavelength of 1000 nm and was absorbed in the center of the strip at initial time t₀ = 0.01 ps. Equation (3.3) was used as initial condition and Eqs. (3.4) and (3.5) were used as boundary conditions.

Figures 4.1 and 4.2 show the electron density C_e(~r, t) and the quasiparticle density C_qp(~r, t), respectively, at several different times after photon absorption. As shown by the simulation, the electron density is continuously spreading throughout the whole strip, with a decreasing concentration at the photon absorption position. This simple behavior is expected since in the diffusion model used, the evolution of C_e(~r, t) follows a diffusion equation without any source term and also it does not depend on C_qp(~r, t). The quasi- particle density is zero at time t₀, it is first being greatly enhanced around the absorption site before it starts to diffuse outwards. Note that the scale of C_e is several orders of mag- nitude smaller than the scale of C_qp, which is reasonable since there is only one excited electron and the total number of exited electron, which is 1, is given by the integral of C_e(~r, t) over the whole strip with the strip area being L × w = 10⁵ nm. Figure 4.3 shows the total number of quasiparticles as a function of time. About 85% of the maximum number of quasiparticles was obtained already after 6 ps.

Next, the results were interpolated onto a Cartesian square grid with a lattice spacing of 0.5 nm in both the x- and the y-direction. The superconducting electron density was calculated based on the result C_qp(~r, t). The non-linear model with pair-breaking effect, Eq. (3.15), was used in the calculation. Supercurrent redistribution was calculated using the iteration method described in section 3.2. Number of iteration steps was set to 1000 before 5 ps and was reduced to 100 steps after that. The bias current was chosen to I_bias = 0.5I_c,v in the x-direction and 0 in the y-direction. Figure 4.4 shows the supercon- ducting electron density at the same chosen moments as for C_e and C_qp. Figures 4.5 and 4.6 are the corresponding x- and y-component of the supercurrent, respectively. Figure

4.7 shows the supercurrent streamline at these moments.

The evolution of n^2D_se(~r, t) is in good match with that of C_qp(~r, t): as shown in Fig.

4.4, the superconducting electron density is successively increasing with time around the strip center, whereas the quasiparticle density is decreasing in that region. This opposite behavior of C_qp and n^2D_se in the center region is expected, since according to the model, the photon absorption should lead to quasiparticle multiplication and thus superconducting electron suppression near the absorption site. Also, since n^2D_se,0 ≈ 0.5 nm⁻², the order of magnitude of n^2D_se should be the same as for C_qp, which it does. Close to the edge regions the superconducting electron density should be less affected and indeed n^2D_se has a higher value and is more homogeneous near the strip edges.

In Fig. 4.5, one can see that the x-component of the supercurrent is suppressed near the center region and enhanced outside that region, which is expected since the supercon- ducting electron density is suppressed in the center so the x-component of supercurrent has to be diverted into the strip side-walks to obey current continuity. In Fig. 4.5d, the supercurrent seems to be going back to a homogeneous state, which is the expected be- havior at long times after the photon absorption. The y-component of the supercurrent is approximately one order of magnitude less than Jx, which is quite reasonable since the bias current we started with was only in the x-direction. Note that J_y has both positive and negative values, implying that the supercurrent goes in different y-directions, some- thing that has to occur in order to preserve current conservation, since the input bias current was zero in the y-direction. Also here the difference tends to smooth out as time increases. Figure 4.7 shows the supercurrent streamlines on the upper half of the strip.

Regions with closer streamlines have higher current densities. Here one can see the effect of current being diverted from the center to the side-walk.

Finally, the potential entry barrier for a single magnetic vortex was calculated. Figure 4.8 compares the integration result using Eq. (3.23) with the analytical solution given by Eq. (3.17), in case of no bias current and uniform bias current throughout the whole strip.

This figure serves as a consistency check since the results should agreed in the absence of any quasiparticles. Figure 4.8a is obtained with a lattice spacing resolution of 0.05 nm. The results are in good agreement. Also, in the absence of bias current, the second term in Eq. (3.17) is zero so the only term left is the logarithmic term that contains a cosine argument, which is an even function. Indeed, the curve of vortex entry barrier in this case is symmetric around the line x = 0. In Fig. 4.8b the bias current was set to Ibias = 0.5Ic,v in the x-direction. Here a good agreement between the analytical and the numerical calculation was achieved already for lattice spacing resolution of 0.2 nm. The vortex core energy is not included in these two plots. However, including it will just intro- duce an overall constant shift of the curves and therefore it will not affect the comparison.

Figure 4.9 shows the vortex entry barrier based on the previous results of supercur- rent and superconducting electron density put into Eq. (3.23). The result was obtained using lattice spacing of 0.5 nm, however, the vortex entry barrier curves are of the same shape and have almost same values as those from simulations with lattice spacing of 0.2 nm. Therefore, the lattice spacing was kept at 0.5 nm in all simulations. The total force exerted on the vortex changes sign at local vortex energy barrier extremum so the

vortex must overcome the barrier hill in order for vortex crossing to take place. There are two potential hills in this case, one near the lower strip edge with the peak around -38 nm and the other one is near the strip center with its peak around 8 nm. The lat- ter one disappears with time, leaving the former as the only potential hill at latest 5.21 ps.

(a) (b)

(c) (d)

Figure 4.1: Electron density C_e(~r, t) for several different times after absorption of a 1000 nm photon in the TaN superconducting strip center. Simulation parameters are listed in Tab. 3.1. Time step was set to dt = 0.1 ps. The evolution of Ce(~r, t) is a simple classic diffusion. Note that the order of magnitude here is 10⁻³ nm⁻², which is reasonable since the strip area is 10⁵ nm² and there is only 1 electron.

(a) (b)

(c) (d)

Figure 4.2: Quasiparticle density C_qp(~r, t) for several different times after absorption of a 1000 nm photon in the TaN superconducting strip center. Simulation parameters are listed in Tab. 3.1. Time step was set to dt = 0.1 ps. Cqp is several orders of magnitude larger than C_e.

Figure 4.3: Number of total quasiparticles in the whole strip as function of time, for absorption of a 1000 nm photon in the TaN strip center. About 85% of the maximum number was obtained after 6 ps.

(a) (b)

(c) (d)

Figure 4.4: Superconducting electron density n^2D_se (~r, t) with calculation based on C_qp(~r, t).

Absorption of a 1000 nm photon in the TaN strip center. ~I_bias = 0.5I_c,vx. The non-linearˆ model with pair-breaking effect was used in the calculation. The value of n^2D_se around the photon absorption site is increasing with time, in accordance with the assumption that local increments of quasiparticles equals local decrements of superconducting electrons and vice versa.

(a) (b)

(c) (d)

Figure 4.5: The x-component of supercurrent after absorption of a 1000 nm photon in the TaN strip center. ~I_bias = 0.5I_c,vx. Current around the absorption site recovers withˆ time in accordance with the increments of super-electron density in this region. The increments of J_x along the side-walks, represented by the dark red regions, is expected due to current conservation. At 10.41 ps, J_x is relaxing back into a homogeneous state.

(a) (b)

(c) (d)

Figure 4.6: The y-component of supercurrent after absorption of a 1000 nm photon in the TaN strip center. ~I_bias = 0.5I_c,vx. The different signs of Jˆ _y indicates different y- directions of the supercurrent, which necessarily has to occur in order to preserve current conservation. J_y relaxes back into a homogeneous state with time.

(a) (b)

(c) (d)

Figure 4.7: Supercurrent streamlines after absorption of a 1000 nm photon in the TaN strip center. ~I_bias = 0.5I_c,vx. Regions with close streamlines correspond to a high super-ˆ current density. The streamlines are "pushed together" from the absorption site into the side-walk. Also, there is a clear tendency of streamlines going back to a homogeneous state especially near the photon absorption site.

(a) (b)

Figure 4.8: Single vortex entry barrier as function of position across the TaN strip, x- coordinate was set to 0. (a) No bias current; (b) Uniform bias current I_bias = 0.5I_c,v in the x-direction. Both case without photon absorption. Numerical calculations were based on Eq. (3.23) excluding the vortex core energy, analytical data points came from Eq. (3.17). The results are in good agreement in both case.

(a) (b)

Figure 4.9: Vortex entry barrier as function of position across the TaN strip, x-coordinate was set to 0. Absorption of a 1000 nm photon at the TaN strip center. I_bias = 0.5I_c,v in the x-direction. (b) is a zoom-in of (a). The vortex has to overcome the potential hills in order to move across the strip. The second hill with its peak around 8 nm is totally disappeared at latest 5.21 ps.

4.1.2 NbN superconducting strip

The simulation was also run with NbN material parameters and the results are shown in Figs. 4.10-4.16.

The evolution of Ce and Cqp share the same characteristic features as in the case with the TaN strip. Here, the total number of quasiparticles reaches about 85% of its maximum value after 6 ps. However, the upper limit of total quasiparticles is about 450 for TaN and only about 300 for NbN. This might not be so strange since the incident photon has the same energy in both case but the superconducting gap is 2.3 meV for NbN and only 1.3 meV for TaN.

The x- and y-component of the supercurrent undergoes same type of change, respec- tively, as in the TaN strip case. Despite this, the evolution of n^2D_se here is very different than the one with TaN strip. In this case, regions to the left and to the right of the center spot have a much higher super-electron density, see the dark red patches in Fig. 4.13.

These two patches move towards the edge and disappear. Also, the value of n^2D_se near the edges of the NbN strip is larger at 10.41 ps than earlier times shown in the same figure.

The vortex entry barrier here spans the same order of magnitude as in the TaN case.

Also here, it has two potential hills at approximately the same positions as those in the TaN case. However, the second potential hill with its peak around 8 nm is much less distinctive than the corresponding one in TaN, and it already disappeared at 2.61 ps.

Also, the shape of the vortex entry barrier curve seems to change less with time for NbN than TaN.

(a) (b)

(c) (d)

Figure 4.10: Electron density C_e(~r, t) for different times after absorption of a 1000 nm photon in the NbN strip center. Simulation parameters are listed in Tab. 3.1. Time step was set to dt = 0.1 ps.

(a) (b)

(c) (d)

Figure 4.11: Quasiparticle density C_qp(~r, t) for different times after absorption of a 1000 nm photon in the NbN strip center. Simulation parameters are listed in Tab. 3.1. Time step was set to dt = 0.1 ps.

Figure 4.12: Number of total quasiparticles in the whole NbN strip, for absorption of a 1000 nm photon in the strip center. About 85% of the maximum number was obtained after 6 ps. The upper limit is only 300, much less than 450 in the TaN case.

(a) (b)

(c) (d)

Figure 4.13: Superconducting electron density n^2D_se(~r, t) with calculation based on C_qp(~r, t). Absorption of a 1000 nm photon in the NbN strip center. ~I_bias = 0.5I_c,vx.ˆ The non-linear model with pair-breaking effect was used in the calculation. Note the dark red patches indicating extra high super-electron density. Also, the value of n^2D_se near strip edges is higher in (d) than in (a) - (c).

(a) (b)

(c) (d)

Figure 4.14: The x-component of supercurrent after absorption of a 1000 nm photon in the NbN strip center. ~I_bias = 0.5I_c,vx. The effect of current suppression around theˆ photon absorption site and thus current enhancement in the side-walks is clearly visible.

(a) (b)

(c) (d)

Figure 4.15: The y-component of supercurrent after absorption of a 1000 nm photon in the NbN strip center. ~I_bias = 0.5I_c,vx. Jˆ _y has both positive and negative values expected from current conservation.

(a) (b)

Figure 4.16: Vortex entry barrier as function of position across the NbN strip, x- coordinate was set to 0. Absorption of a 1000 nm photon at the NbN strip center.

I_bias = 0.5I_c,v in the x-direction. (b) is a zoom-in of (a). The second potential hill is much less distinctive in shape and disappears earlier than the corresponding one in the TaN strip case. Also, the shape of the potential entry barrier changes less with time here compared with TaN.

4.2 Photon absorption near strip edge

This section shows the simulation results in case of photon absorption near the lower strip edge: x = 0, y = −w/2 + 1.5ξ. All other simulation parameters are kept the same as in section 4.1.

4.2.1 TaN superconducting strip

Simulation results are shown in Figs. 4.17-4.23.

In this case the quasiparticle density and the superconducting electron density seems to build up at the lower edge after diffusion from the absorption site. This "building up" process dominates over the outward diffusion at the lower edge between 2.61 ps and 5.21 ps. This probably occurs due to the Neumann boundary condition. Again, the total number of quasiparticles reached about 85% of its maximum value after only 6 ps. The upper limit in this case is about 375, which is a bit less than 450, which is the upper limit of total quasiparticles in case of photon absorption at TaN strip center. The reason for this is probably also the Neumann boundary condition.

The x-component of the supercurrent again shows the characteristic behavior: sup- pression around the photon absorption site, which in this case is at y ≈ −42 nm, and enhancement on the side-walk. The y-component of the supercurrent takes both positive and negative values as usual. However, unlike photon absorption at the strip center, there is no clear sign of supercurrent relaxing back into a homogeneous state at 10.41 ps.

The main difference of the vortex entry barrier between absorption near lower edge and absorption in center is the disappearance of the second potential hill in the former case, compare Fig. 4.23a with Fig. 4.9a. Another thing to notice is that for photon ab- sorption at strip center, the value of the potential decreases with time on the lower half of the strip and increases with time on the upper half of the strip, whereas in this case the overall potential barrier increases with time. Also, the shape of the potential curve does not change as much here as the one with photon absorption at center.

(a) (b)

(c) (d)

Figure 4.17: C_e(~r, t) for different times after absorption of a 1000 nm photon at 1.5 coherence length away from the TaN lower strip edge.

(a) (b)

(c) (d)

Figure 4.18: C_qp(~r, t) for different times after absorption of a 1000 nm photon at 1.5 coherence length away from the TaN lower strip edge. C_qp(~r, t) is piling up at the lower strip edge in (b) - (c). This is probably due to the imposed Neumann boundary condition, which says that particles can not pass through the strip edge.

Figure 4.19: Number of total quasiparticles in the whole TaN strip as function of time.

Absorption of a 1000 nm photon near lower strip edge. More than 85% of the maximum quasiparticle number was reached after 6 ps. The upper limit in this case is about 375, whereas for photon absorption at TaN strip center the upper limit is about 450.

(a) (b)

(c) (d)

Figure 4.20: n^2D_se (~r, t) with calculation based on C_qp(~r, t). Absorption of a 1000 nm photon at 1.5 coherence length away from the TaN lower strip edge. ~I_bias = 0.5I_c,vx. Theˆ non-linear model with pair-breaking effect was used. n^2D_se (~r, t) is piling up at the lower edge in (b) - (c), corresponding to the same effect caused by the boundary condition on C_qp(~r, t) shown in Fig. 4.18.

(a) (b)

(c) (d)

Figure 4.21: The x-component of supercurrent at several different times after absorption of a 1000 nm photon near the TaN lower strip edge. ~I_bias = 0.5I_c,vx. In this case, Jˆ _x does not seem to relax back into a homogeneous value 10.41 ps after photon absorption.

(a) (b)

(c) (d)

Figure 4.22: The y-component of supercurrent at several different times after absorption of a 1000 nm photon near the TaN lower strip edge. ~I_bias = 0.5I_c,vx. Compared to theˆ case with photon absorption at the TaN strip center, in this case there is no clear sign that J_y is relaxing back into a homogeneous value throughout the strip.

(a) (b)

Figure 4.23: Vortex entry barrier as function of position across the strip, x-coordinate was set to 0. I_bias = 0.5I_c,v in the x-direction. Absorption of a 1000 nm photon 1.5 coherence length away from the TaN lower strip edge. (b) is a zoom-in of (a). In this case, there is only one potential hill and overall the shape of the curve stays approximately the same.

Overall, the value of potential entry barrier is shifting slightly upwards with time.

4.2.2 NbN superconducting strip

The simulation results are shown in Figs. 4.24-4.30.

The evolution of C_e and C_qp show similar features as those for the TaN strip. More than 85% of the total quasiparticle upper limit was again reached after 6 ps. The value of quasiparticle upper limit in this case is about 200, which is approximately 30% less than the upper limit for photon absorption at NbN strip center and 45% less than the value in the case of photon absorption at TaN edge.

Simulations results for n^2D_se (~r, t) shows that regions to the left and to the right of the photon absorption spot have extra high superconducting electron density, just as in the case of photon striking at NbN strip center shown in Fig. 4.13. The x-component of the supercurrent is again suppressed around the photon absorption spot and enhanced in the side-walk. The relaxation of J_x back into a homogeneous state throughout the whole strip seems to happen faster here than in the case of photon striking near TaN strip edge.

However, the relaxation here is not as fast as for photon absorption at NbN strip center.

Regarding to J_y, one can not see relaxation happening before 10.41 ps.

Like the case with photon absorption near TaN strip edge, the potential entry bar- rier has only one hill here, located near the lower strip edge. Also, overall the value of the potential is slightly increased with time, and keeps approximately the same shape.

However, the height of the potential barrier is much higher in this case, with the highest peak at 1.01 ps of the four times shown, whereas the highest peak for TaN is at 10.41 ps.

(a) (b)

(c) (d)

Figure 4.24: C_e(~r, t) for different times after absorption of a 1000 nm photon at 1.5 coherence length away from the NbN lower strip edge.

(a) (b)

(c) (d)

Figure 4.25: C_qp(~r, t) for different times after absorption of a 1000 nm photon at 1.5 coherence length away from the NbN lower strip edge.

Figure 4.26: Number of total quasiparticles in the whole NbN strip as function of time.

Absorption of a 1000 nm photon near lower strip edge. More than 85% of maximum number was reached after 6 ps. The upper limit here is about 200, which is approximately 100 less than photon absorption at NbN strip center, and approximately 175 less than photon absorption near TaN strip edge.

(a) (b)

(c) (d)

Figure 4.27: n^2D_se (~r, t) with calculation based on C_qp(~r, t). Absorption of a 1000 nm photon at 1.5 coherence length away from the NbN lower strip edge. ~I_bias = 0.5I_c,vx.ˆ The non-linear model with pair-breaking effect was used. As in the case with photon absorption at NbN strip center, there are regions outside the photon absorption spot with extra high super-electron density.

(a) (b)

(c) (d)

Figure 4.28: The x-component of supercurrent after absorption of a 1000 nm photon near the NbN lower strip edge. ~I_bias = 0.5I_c,vx. The characteristic features of supercur-ˆ rent suppression and supercurrent enhancement in different regions are clearly visible.

Relaxation of J_x back into a homogeneous state takes longer time here compared with the corresponding process for photon absorption at NbN strip center shown in Fig. 4.14, but it happens sooner compared to photon absorption near TaN strip edge shown in Fig.

4.21.

(a) (b)

(c) (d)

Figure 4.29: The y-component of supercurrent after absorption of a 1000 nm photon near the NbN lower strip edge. ~I_bias = 0.5I_c,vx. As in the case of photon absorption near TaNˆ strip edge, Jy takes both positive and negative values just as expected, but there is no clear sign of it relaxing back into a homogeneous value throughout the strip.

(a) (b)

Figure 4.30: Vortex entry barrier as function of position across the strip, x-coordinate was set to 0. I_bias = 0.5I_c,v in the x-direction. Absorption of a 1000 nm photon 1.5 coherence length away from the TaN lower strip edge. (b) is a zoom-in of (a). As in the case of photon striking near TaN strip edge, there is only one potential hill here, and overall the potential curve is moving slightly upwards with time without much change in its shape.

Chapter 5 Summary

The photon detection process in SNSPD can be divided into three steps: photon ab- sorption, normal domain generation, and electrical signal output [4]. The first and the last step of this process have been improved to significantly increase the overall pho- ton detection efficiency. However, the mechanism responsible for the second step is not fully understood [4]. Therefore, the aim of this thesis is to gain more understanding of the normal-conducting domain generation process. Following Ref. [5], this process was further divided into three steps: quasiparticle multiplication and diffusion, supercurrent redistribution, and triggering of voltage pulse. The theory behind all these three steps is explained in Ch. 3. During the investigations, a more suitable approach for calculat- ing supercurrent redistribution was found: instead of only requiring current continuity, we also added the condition of supercurrent kinetic energy minimization. Besides this, the contribution of vortex core energy was added to the vortex entry barrier for single vortex crossing through the superconducting strip. Motivations behind these modifica- tions can be found in the corresponding sections above. Simulations were done for two cases: photon absorption at the center of the superconducting strip and absorption at 1.5 coherence length away from the strip lower edge. In each case, the simulation was done for both TaN and NbN superconducting strip. Simulation results are shown in Ch. 4.

It was found that in all cases considered, about or more than 85% of the maximum number of quasiparticles was reached after 6 ps. For each material, the maximum number of quasiparticles is higher for photon absorption at the strip center than absorption near the strip edge. The relaxation of supercurrent back into a homogeneous state throughout the strip takes shorter time in case of photon absorption at the strip center.

One interesting phenomenon that appeared during the simulations is the formation of regions with extra high super-electron density to the left and to the right of the photon absorption site for the case of a NbN superconducting strip, both in the case of photon absorption at strip center, see Fig. 4.13, and in the case of photon absorption near strip edge, see Fig. 4.27. This phenomenon did not occur in simulations of a TaN supercon- ducting strip.

It was also found that for both materials, the vortex entry barrier in case of photon absorption at the strip center had two local maxima, one located near the strip lower edge and the other near the strip center. The latter maximum disappeared as time increased.

For photon absorption at 1.5 coherence length away from the strip edge, the vortex en- try barrier only showed one local maximum near the lower edge for both materials. In this case, the overall shape of the potential curve stayed approximately the same and only shifted slightly upwards as time increased. However, in order to obtain more ac- curate results for the potential entry barrier a more refined model needs to be considered.

One major refinement would be to implement the motion of vortex crossing through the strip width while taking into account the effect of current crowding outside the vortex core. One complication concerning this refinement is that to simulate vortex motion one needs a microscopic theory to accurately model the vortex core region. The other com- plication is that since supercurrent redistribution, which is calculated using an iterative method, has to be performed at vortex position through the strip width, the simulation will be computationally expensive. However, this might be necessary for improving the photon detection model and can be considered as a next step to the present investigation.

There are several simplifications made in the model. As pointed out in article [5], both the presence of quasiparticles and bias current lead to a certain reduction of the superconducting gap ∆. Also, it is not entirely correct to assume a constant time scale τqp for the quasiparticle multiplication process. In addition, the diffusion coefficient may be a function of the excitation energy of the electron and thus not constant in time. Also, in the simulations done, thermal fluctuations were ignored.

In conclusion, this project suggests some improvements of the current model for pho- ton detection in SNSPD and shows simulation results based on the modifications made.

One major refinement is suggested as a next step of investigation to further improve the model. Hopefully, this project takes us one step closer towards the goal of fully under- standing of the photon detection mechanism of SNSPD.