Grand canonical partition function of a serial metallic island system

A circuit diagram of the finite serial island system is depicted in Fig.1. Microscopically, this system can be theoretically modeled by the Hamiltonian [1]

$\displaystyle H=H_{B}+H_{T}+H_{C}\,,$

(1)

where the term $H_{B}$ describes the conduction electrons of the leads and the islands

$\displaystyle H_{B}=\sum_{Jk\sigma}\epsilon_{Jk\sigma}c_{Jk\sigma}^{{\dagger}}c_{Jk\sigma}+\sum_{Ik\sigma}\epsilon_{Ik\sigma}d_{Ik\sigma}^{{\dagger}}d_{Ik\sigma}\,,$

(2)

respectively. Here $\epsilon_{Jk\sigma}$ is the energy of an electron with longitudinal wave vector $k$ in channel $\sigma$ of lead $J$, where $J\in\{S,D\}$. The channel index $\sigma$ includes the transversal and spin quantum numbers. $c_{Jk\sigma}$ and $c_{Jk\sigma}^{{\dagger}}$ are the corresponding electron annihilation and creation operators, respectively. Likewise, $\epsilon_{Ik\sigma}$ denotes the energy of an electron with the longitudinal wave vector $k$ in channel $\sigma$ of island $I$, where $I\in\{1,2,...,N\}$, and $d_{Ik\sigma}$ and $d_{Ik\sigma}^{{\dagger}}$ are the associated electron annihilation and creation operators, respectively.

The second term in Eq.(1) describes the tunneling of electrons across the finite tunneling junctions of the circuit as shown in Fig.1,

	$\displaystyle H_{T}=\sum_{kq\sigma}$	$\displaystyle\Big{[}$	$\displaystyle d_{1k\sigma}^{{\dagger}}t_{1Skq\sigma}{\rm e}^{-i\varphi_{1}}c_{Sq\sigma}+\sum_{I=2}^{N}d_{Ik\sigma}^{{\dagger}}t_{II-1kq\sigma}{\rm e}^{-i(\varphi_{I}-\varphi_{I-1})}c_{I-1q\sigma}$
			$\displaystyle+c_{Dk\sigma}^{{\dagger}}t_{DNkq\sigma}{\rm e}^{i\varphi_{N}}d_{Nq\sigma}+\hbox{H.c.}\Big{]},$

where $t_{1Skq\sigma}$, $t_{II-1kq\sigma}$, and $t_{DNkq\sigma}$ denote electron tunneling amplitude of the tunneling junctions labeled by $1S,II-1$, and $DN$, respectively. For example, $t_{1Skq\sigma}$ is the amplitude for an electron in state $|q\sigma\rangle$ on lead $S$ to tunnel onto island $1$ with final state $|k\sigma\rangle$. In addition, we have introduced the phase operators $\varphi_{I}$ conjugate to the number operators $n_{I}$ of excess charges on the islands $I$, i.e. $[n_{I},\varphi_{I}]=i$. Accordingly, the charge shift operator ${\rm e}^{-i\varphi_{I}}$ add one charge to the island $I$. Here and in the following, we have defined $\varphi_{S}=\varphi_{D}=0$ since $V_{S}=V_{D}=0$.

The last term in Eq.(1) is Coulomb charging energy expressed as[1, 24]

$\displaystyle H_{C}=\sum\limits_{I=1}^{N}{{E_{II}}{({n_{I}}-{n_{0I}})}^{2}}+2\sum\limits_{I<I^{\prime}}^{N}{{E_{II^{\prime}}}({n_{I}}-{n_{0I}})({n_{I^{\prime}}}-{n_{0I^{\prime}}})},$

(4)

where $n_{I}$ is the (excess) electron number in the island $I$. For the island $I$ the continuous charge induced by the gate voltage $I$ is denoted by $n_{0I}$ expressed by $n_{0I}=C_{gI}V_{gI}/e$. The variable $n_{0I}$ is non-integer and just a parameter of the external voltage. The coefficients $E_{II}$ and $E_{II^{\prime}}$ are the matrix elements of the matrix $\mathbf{E}_{CN}$ defined by

$\displaystyle\mathbf{E}_{CN}=\frac{{{e}^{2}}}{2}\mathbf{C}_{N}^{-1},$

(5)

where $\mathbf{C}_{N}^{-1}$ is the inverse matrix of $\mathbf{C}_{N}$

$\displaystyle\mathbf{C}_{N}=\begin{pmatrix}C_{11}&C_{12}&\cdots&C_{1N}\\ C_{21}&C_{22}&\cdots&C_{2N}\\ \vdots&\vdots&\vdots&\vdots\\ C_{N1}&\cdots&\cdots&C_{NN}\end{pmatrix},$

(6)

and the matrix elements of $C_{II^{\prime}}$ are defined by the conditions:

$\displaystyle C_{II^{\prime}}=\begin{cases}0,&\text{if}\,\,\,\,\,|I-I^{\prime}|>1\\ C_{\sum\nolimits I},&\text{if}\,\,\,\,\,I=I^{\prime}\\ -C_{I^{\prime}},&\text{if}\,\,\,\,\,|I-I^{\prime}|=1,\end{cases}$

(7)

where $C_{\sum\nolimits I}$ is the sum of all capacitors connected with the island $I$, i.e., $C_{\sum I}=C_{I}+C_{I-1}+C_{gI}$. According to the conditions in Eq.(7), one obtains $C_{II^{\prime}}=C_{I^{\prime}I}$.

The goal of this section is a path integral representation of the grand partition function

$\displaystyle Z={\rm tr}\,\{{\rm e}^{-\beta(\hat{H}-\mu\hat{N})}\}$

(8)

of the island system suitable for the evaluation, where $\mu$ denotes the chemical potential of the system and $\hat{N}$ is the particle number operator. The Hilbert space of the Hamiltonian in Eqs.(1)–(4) is the product of the space spanned by the state $\ket{n}$, or equivalently the phase state $\ket{\varphi}$, and the Fock space of the quasi–particles. In order to trace over the quasi–particle excitation, we introduce Grassmann numbers $\zeta_{Jk\sigma}$ and $\theta_{Ik\sigma}$ corresponding to the electron annihilators $c_{Jk\sigma}$ and $d_{Ik\sigma}$, respectively. For notation convenience, the Grassmann numbers for electronic lead states are combined to a single vector $\vec{\zeta}=(\zeta_{Sk\sigma},\zeta_{Dq\sigma})^{T}$. Likewise, $\vec{\theta}=(\theta_{1q\sigma},\theta_{2q\sigma},\ldots,\theta_{Nq\sigma})^{T}$ combines all Grassmann numbers for electronic island states. To trace over the island charges, we work in the phase representation and introduce the vector $\vec{\varphi}=(\varphi_{1},\ldots,\varphi_{N})^{T}$. The partition function in Eq.(8) may then be written as

$\displaystyle Z=\int D\mu(\vec{\varphi})\int D\mu(\vec{\zeta})\int D\mu(\vec{\theta}){\rm e}^{-\vec{\zeta}^{\ast}\cdot\vec{\zeta}-\vec{\theta}^{\ast}\cdot\vec{\theta}}\langle\vec{\varphi},-\vec{\zeta},-\vec{\theta}|\,{\rm e}^{-\beta H}\,|\vec{\varphi},\vec{\zeta},\vec{\theta}\rangle\,,$

(9)

where we have dropped the term $(\mu\hat{N})$ in the partition function. This term is further discussed this term in Section 2.4. Here $|\vec{\varphi}\rangle$ stands for the product of phase states $|\varphi_{I}\rangle$. The shorthand notations $|\vec{\zeta}\rangle$ and $|\vec{\theta}\rangle$ stand for the product of fermion coherent states $|\zeta_{Jk\sigma}\rangle$ and $|\theta_{Ik\sigma}\rangle$, respectively. The shorthand integration

$\displaystyle\int D\mu(\vec{\zeta})\equiv\int\prod_{Jk\sigma}d\zeta_{Jk\sigma}^{\ast}\,d\zeta_{Jk\sigma},$

(10)

and

$\displaystyle\int D\mu(\vec{\theta})\equiv\int\prod_{Iq\sigma}d\theta_{Iq\sigma}^{\ast}\,d\theta_{Iq\sigma}$

(11)

are over all Grassmann numbers related with the lead and island states, respectively. The short hand integration

$\displaystyle\int D\mu(\vec{\varphi})\,\equiv\,\left(\frac{1}{2\pi}\right)^{2}\int\prod_{I}\,d\varphi_{I}$

(12)

integrates over the phases $\varphi_{I}$ that are $2\pi$-periodic. Further, we have introduced the notations

$\displaystyle\vec{\zeta}^{\ast}\cdot\vec{\zeta}=\sum_{Jk\sigma}\zeta^{\ast}_{Jk\sigma}\zeta_{Jk\sigma},$

(13)

and

$\displaystyle\vec{\theta}^{\ast}\cdot\vec{\theta}=\sum_{Iq\sigma}\theta^{\ast}_{Iq\sigma}\theta_{Iq\sigma}\,.$

(14)

The derivation of the path integral expression can be done as usual by multiple insertions of the closure relation [21]

$\displaystyle 1=\int D\mu(\vec{\varphi})\int D\mu(\vec{\zeta})\int D\mu(\vec{\theta}){\rm e}^{-\vec{\zeta}^{\ast}\cdot\vec{\zeta}-\vec{\theta}^{\ast}\cdot\vec{\theta}}\,\ket{\vec{\varphi},\vec{\zeta},\vec{\theta}}\bra{\vec{\varphi},\vec{\zeta},\vec{\theta}},$

(15)

in the product space. The imaginary–time step is introduced as $\Delta=\beta/P$, where we have used $\hbar=1$ through this paper, and $P$ denotes the Trotter number. For each imaginary time segment $\Delta_{j}=\tau_{j}-\tau_{j-1},$ one can calculate the short-time propagator

	$\displaystyle\langle\vec{\varphi}_{j},\vec{\zeta}_{j},\vec{\theta}_{j}|\,{\rm e}^{-\Delta_{j}H}\,|\vec{\varphi}_{j-1},\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle$	$\displaystyle=$	$\displaystyle\langle\vec{\zeta}_{j},\vec{\theta}_{j}|\,{\rm e}^{-\Delta_{j}\left[H_{B}+H_{T}(\vec{\varphi}_{j-1})\right]}\,|\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle$		(16)
			$\displaystyle\times\langle\vec{\varphi}_{j}|{\rm e}^{-\Delta_{j}H_{C}}|\vec{\varphi}_{j-1}\rangle,$

where $\vec{\varphi}_{j}=(\varphi_{1}(\tau_{j}),\varphi_{2}(\tau_{j}),\ldots,\varphi_{N}(\tau_{j}))^{T}$ and $\vec{\varphi}_{j-1}=(\varphi_{1}(\tau_{j-1}),\varphi_{2}(\tau_{j-1}),\ldots,\varphi_{N}(\tau_{j-1}))^{T}$ and $\vec{\zeta}_{j}$ and $\vec{\theta}_{j}$ are the Grassmann variables at time $\tau_{j}$. Further, $H_{T}(\vec{\varphi}_{j})$ refers to the Hamiltonian in Eq.(2.1) with the phase operators $\varphi_{I}$ replaced by $\varphi_{I}(\tau_{j})$. Making use of the form in Eq.(16) of the short time propagator in the Trotter break-up for the partition function, we obtain

	$\displaystyle Z$	$\displaystyle=$	$\displaystyle\int\limits_{0}^{2\pi}\prod_{I=1}^{N}\Big{[}d\varphi_{I,P}\ldots d\varphi_{I,0}\,\delta(\varphi_{I,P}-\varphi_{I,0})\Big{]}$
			$\displaystyle\times\langle\vec{\varphi}_{P}|\,{\rm e}^{-\Delta_{P}H_{C}}\,|\vec{\varphi}_{P-1}\rangle\ldots\langle\vec{\varphi}_{1}|\,{\rm e}^{-\Delta_{1}H_{C}}\,|\vec{\varphi}_{0}\rangle\,Z_{BT}[\vec{\varphi}],$

where $Z_{BT}[\vec{\varphi}]$ is the trace over the Grassmann fields discussed in more detail in Section 2.4.

The representation of the partition function in Eq.(2.2) can be evaluated by inserting the closure relation

$\displaystyle 1=\int\limits_{0}^{2\pi}{d\varphi\left|\varphi\right\rangle}\left\langle\varphi\right|=\sum\limits_{n=-\infty}^{\infty}{\left|n\right\rangle}\left\langle n\right|,$

(18)

where $n$ is an integer number into the short time propagator in the last term in Eq.(16). We have

	$\displaystyle\langle\vec{\varphi}_{j}|{\rm e}^{-\Delta_{j}H_{C}}|\vec{\varphi}_{j-1}\rangle$	$\displaystyle=$	$\displaystyle\sum_{\vec{n}}\langle\varphi_{1,j}|n_{1}\rangle\,\ldots\langle\varphi_{N,j}|n_{N}\rangle$
			$\displaystyle\times{\rm e}^{-\Delta_{j}H_{C}(n_{1},\ldots,n_{N})}\langle n_{1}|\varphi_{1,j-1}\rangle\ldots\langle n_{N}|\varphi_{N,j-1}\rangle$
		$\displaystyle=$	$\displaystyle\frac{1}{({2\pi})^{N}}\sum_{\vec{n}}{\rm e}^{-\Delta_{j}H_{C}(n_{1},\ldots,n_{N})}{\rm e}^{-i\sum_{I}n_{I}(\varphi_{I,j}-\varphi_{I,j-1})},$

where we have defined

$\displaystyle\sum_{\vec{n}}\equiv\sum_{n_{1}=-\infty}^{\infty}...\sum_{n_{N}=-\infty}^{\infty}.$

(20)

$H_{C}(n_{1},\ldots,n_{N})$ is the Coulomb Hamiltonian in Eq.(4) for given eigenvalues of the electron number operators $n_{I}$ and $\left\langle{{{n_{I}}}}\mathrel{\left|{\vphantom{{{n_{I}}}{{\varphi_{I,j}}}}}\right.\kern-1.2pt}{{{\varphi_{I,j}}}}\right\rangle={\left({2\pi}\right)^{-1/2}}\exp\left({i{n_{I}}{\varphi_{I,j}}}\right)$. By means of the Poisson re-summation formula [21],

$\displaystyle\sum_{n=-\infty}^{\infty}f(n)=\sum_{k=-\infty}^{\infty}\ \int\limits_{-\infty}^{\infty}\!\!dn\,{\rm e}^{-2\pi ikn}f(n)\,,$

(21)

and the matrix in Eq.(5) corresponding with $H_{C}$, the short time propagator can be transformed to read

	$\displaystyle\langle\vec{\varphi}_{j}|{\rm e}^{-\Delta_{j}H_{C}}|\vec{\varphi}_{j-1}\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{(2\pi)^{N}}\sum_{\vec{k}_{j}}\int\limits_{-\infty}^{\infty}d\tilde{n}_{1}\ldots\int\limits_{-\infty}^{\infty}d\tilde{n}_{N}$
			$\displaystyle\times\exp\left[{-i\Delta_{j}\sum_{I=1}^{N}\tilde{n}_{I}\Delta\varphi_{I}}-{\Delta_{j}\vec{\tilde{n}}_{I}^{T}\mathbf{E}_{CN}\vec{\tilde{n}}_{I}}\right]\,,$

where we have defined

$\displaystyle\sum_{\vec{k}_{j}}\equiv\sum_{k_{1,j}=-\infty}^{\infty}...\sum_{k_{N,j}=-\infty}^{\infty},$

(23)

$\vec{\tilde{n}}_{I}={\left({{{\tilde{n}}_{1}},...,{{\tilde{n}}_{N}}}\right)^{T}}$ with $\tilde{n}_{I}\equiv(n_{I}-n_{I0})$, and $\Delta{\varphi_{I}}=\left({{\varphi_{I,j}}-{\varphi_{I,j-1}}+2\pi{k_{I,j}}}\right)/{\Delta_{j}}$. Evaluating the Gaussian integrals for matrix form [21]

	$\displaystyle\int\limits_{-\infty}^{\infty}{d{n_{1}}...d{n_{N}}}\exp[{-\lambda({{\vec{n}^{T}}{\mathbf{E}_{CN}}\vec{n}+{\vec{n}^{T}}\vec{J}})}]$	$\displaystyle=$	$\displaystyle{({\frac{\pi}{\lambda}})^{N/2}}[{\det{\mathbf{E}_{CN}}]^{-\frac{1}{2}}}$
			$\displaystyle\times\exp[{({\frac{\lambda}{4}}){\vec{J}}^{T}{\mathbf{E}_{CN}^{-1}}\vec{J}}],$

where $\vec{n}={\left({{{\tilde{n}}_{1}},...,{{\tilde{n}}_{N}}}\right)^{T}}$, $\vec{J}={\left({\Delta{\varphi_{1}},...,\Delta{\varphi_{N}}}\right)^{T}}$, and $\lambda$ is a constant, one obtains

$\displaystyle\langle\vec{\varphi}_{j}|{\rm e}^{-\Delta_{j}H_{C}}|\vec{\varphi}_{j-1}\rangle={\mathcal{N}_{j}}\sum_{\vec{k}_{j}}{\exp\left[{\frac{{-{\Delta_{j}}}}{4}\Delta\vec{\varphi}_{j}^{T}\mathbf{E}_{CN}^{-1}\Delta{\vec{\varphi}_{j}}-i{\Delta_{j}}\vec{\bm{n}}_{g}^{T}\Delta{\vec{\varphi}_{j}}}\right]}$

(25)

where $\mathcal{N}_{j}$ is the normalization constant for the time segment $j$;

$\displaystyle{\mathcal{N}_{j}}=\frac{1}{{{{\left({2\pi}\right)}^{N}}}}{\left({\frac{\pi}{{{\Delta_{j}}}}}\right)^{N/2}}{\left[{\det{\mathbf{E}_{CN}}}\right]^{-\frac{1}{2}}},$

(26)

that can be incorporated into the path integral measure. We now use the freedom to relabel the summations over the winding numbers $k_{I}$ and to transform the integrals over $\varphi_{I}$. Instead of summations over $k_{I,j},\ldots,k_{I,P}$, we sum over

$\displaystyle k^{\prime}_{I,n}=\sum_{j=1}^{n}k_{I,j}\quad\mbox{for}\quad n=1,\ldots,P\,,$

(27)

with the consequence that

$\displaystyle\Delta\varphi_{I,j}=\frac{\varphi_{I,j}+2\pi k^{\prime}_{I,j}-\varphi_{I,j-1}-2\pi k^{\prime}_{I,j-1}}{\Delta_{j}}.$

(28)

Using $\varphi^{\prime}_{I,j}=\varphi_{I,j}+2\pi k^{\prime}_{I,j}$ and dropping the primes as a convenient integration variable, we can rewrite the partition function in Eq.(2.2) as

	$\displaystyle Z$	$\displaystyle=$	$\displaystyle\mathcal{N}\prod\limits_{I=1}^{N}\Bigg{[}\sum_{k_{I,j}=-\infty}^{\infty}\int\limits_{2\pi k_{I,P}}^{2\pi(k_{I,P}+1)}d\varphi_{I,P}\ldots\int\limits_{2\pi k_{I,1}}^{2\pi(k_{I,1}+1)}d\varphi_{I,1}\int\limits_{0}^{2\pi}d\varphi_{I,0}$
			$\displaystyle\times\delta(\varphi_{I,P}-\varphi_{I,0}-2\pi k_{I,P})\Bigg{]}$
			$\displaystyle\times\exp[\sum_{j=1}^{P}\Delta_{j}(\frac{1}{4}\Delta\vec{\varphi}^{T}_{j}\mathbf{E}_{CN}^{-1}\Delta\vec{\varphi}_{j}+i\vec{n}^{T}_{g}\cdot\Delta\vec{\varphi}_{j})]Z_{BT}[\vec{\varphi}],$

where we have introduced the vectors

$\displaystyle\Delta\vec{\varphi}_{j}=\left(\frac{\varphi_{1,j}-\varphi_{1,j-1}}{\Delta_{j}},\frac{\varphi_{2,j}-\varphi_{2,j-1}}{\Delta_{j}},...,\frac{\varphi_{N,j}-\varphi_{N,j-1}}{\Delta_{j}}\right)^{T},$

(30)

and $\vec{n}_{g}=\left(n_{01},n_{02},...n_{0N}\right)^{T}$. The normalization factor $\mathcal{N}\equiv\prod_{j=1}^{P-1}\mathcal{N}_{j}$.

With the exception of $k_{P}$, the sums over $k_{j}$ can be incorporated in the integrals over $\varphi_{I,j}$. This simplifies the expression in Eq.(2.3) further, and we get in the continuum limit, i.e., ($P\rightarrow\infty,\ \Delta_{j}\rightarrow 0$) which is a path integral over all imaginary-time paths $\vec{\varphi}(\tau)=(\varphi_{1}(\tau),\varphi_{2}(\tau),...,\varphi_{N}(\tau))^{T}$ in the time interval $(0,\beta)$ with arbitrary winding numbers

$\displaystyle Z$

$\displaystyle=$

$\displaystyle\mathcal{N}\prod\limits_{I=1}^{N}\bigg{[}\sum_{k_{I}=-\infty}^{\infty}\int\limits_{\varphi_{I}(0)}^{\varphi_{I}(\beta)+2\pi k_{I}}D[\varphi_{I}(\tau)]\bigg{]}{\rm e}^{-S_{C}[\vec{\varphi}(\tau)]}\,Z_{BT}[\vec{\varphi}]\,,$

(31)

and the Coulomb action of the serial island system is expressed by

$\displaystyle S_{C}[\vec{\varphi}(\tau)]=\int_{0}^{\beta}d\tau\left[\frac{1}{4}\dot{\vec{\varphi}}^{T}\mathbf{E}_{CN}^{-1}\,\dot{\vec{\varphi}}+i(\vec{n}^{T}_{g}\cdot\dot{\vec{\varphi}})\right]\,,$

(32)

where $\dot{\vec{\varphi}}(\tau_{j})=\left(\dot{\varphi}_{1}(\tau_{j}\right),\dot{\varphi}_{2}(\tau_{j}),...,\dot{\varphi}_{N}(\tau_{j}))^{T}$ is the continuum limit of $\Delta\vec{\varphi}_{j}$ for $\Delta_{j}\rightarrow 0$ and $\mathbf{E}_{CN}^{-1}$ is the inverse matrix of the matrix $\mathbf{E}_{CN}$ in Eq.(5).

In this section, we have to evaluate the partition function in Eq.(31) wherein $Z_{BT}$ involves the integration over the Grassmann numbers. In the Trotter break-up discussed in the previous section, the imaginary-time segment $\Delta_{j}$ gives rise to the short time propagator in Eq.(16) with the fermions short time propagator

$\displaystyle\langle\vec{\zeta}_{j},\vec{\theta}_{j}|\,{\rm e}^{-\Delta_{j}\left[H_{B}+H_{T}(\vec{\varphi}_{j-1})\right]}\,|\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle.$

(33)

Within this short time propagator, the electron annihilation (creation) operators in the Hamiltonians $H_{B}$ in Eq.(2) and $H_{T}$ in Eq.(2.1) can be replaced for small $\Delta_{j}$ by the corresponding (conjugate) Grassmann numbers. Accordingly, one finds

	$\displaystyle\langle\vec{\zeta}_{j},\vec{\theta}_{j}|\,{\rm e}^{-\Delta_{j}\left[H_{B}+H_{T}(\vec{\varphi}_{j-1})\right]}\,|\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle$	$\displaystyle=\langle\vec{\zeta}_{j},\vec{\theta}_{j}|\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle$		(34)
		$\displaystyle\times\exp\{-\Delta_{j}\Big{[}H_{B,j,j-1}+H_{T,j,j-1}(\vec{\varphi}_{j-1})\Big{]}\bigg{\}},$

where

$\displaystyle H_{B,j,j-1}=\sum_{Jk\sigma}\epsilon_{Jk\sigma}\zeta^{*}_{Jk\sigma,j}\zeta_{Jk\sigma,j-1}+\sum_{Ik\sigma}\epsilon_{Ik\sigma}\theta^{*}_{Ik\sigma,j}\theta_{Ik\sigma,j-1}\,,$

(35)

and

	$\displaystyle H_{T,j,j-1}(\vec{\varphi}_{j-1})$	$\displaystyle=$	$\displaystyle\sum_{kq\sigma}\Big{[}\theta_{1k\sigma,j}^{*}t_{1Skq\sigma}\,{\rm e}^{-i\varphi_{1,j-1}}\zeta_{Sq\sigma,j-1}$
			$\displaystyle+\sum_{I=2}^{N}\theta_{Ik\sigma,j}^{*}t_{II-1kq\sigma}\,{\rm e}^{-i(\varphi_{I,j-1}-\varphi_{I-1,j-1})}\theta_{I-1q\sigma,j-1}$
			$\displaystyle+\zeta_{Dk\sigma,j}^{*}t_{DNkq\sigma}\,{\rm e}^{i\varphi_{N,j-1}}\theta_{Nq\sigma,j-1}+\hbox{H.c.}\Big{]}\,.$

Combining the exponential factor in Eq.(34) with the factor $\exp\left(-\vec{\zeta}^{*}_{j}\cdot\vec{\zeta}_{j}-\vec{\theta}_{j}^{*}\cdot\vec{\theta}_{j}\right)$ stemming from the closure relation in Eq.(15), and the scalar product
$\langle\vec{\zeta}_{j},\vec{\theta}_{j}|\vec{\zeta}_{j-1},\vec{\theta}_{j-1}\rangle=\exp\left(\vec{\zeta}_{j}^{*}.\vec{\zeta}_{j-1}+\vec{\theta}_{j}^{*}.\vec{\theta}_{j-1}\right)$, we obtain for each time segment $\Delta_{j}$ a factor

$$\exp\{-\Delta_{j}\bigg{[}\vec{\zeta}_{j}^{*}\cdot\frac{\vec{\zeta}_{j}-\vec{\zeta}_{j-1}}{\Delta_{j}}+\vec{\theta}^{*}_{j}\cdot\frac{\vec{\theta}_{j}-\vec{\theta}_{j-1}}{\Delta_{j}}+H_{B,j,j-1}+H_{T,j,j-1}(\vec{\varphi}_{j-1})\bigg{]}\bigg{\}}\,.$$

(37)

In the continuum limit, these factors and the integration over the Grassmann numbers $\vec{\zeta}_{j}$ and $\vec{\theta}_{j}$ are combined to the path integral

$\displaystyle Z_{BT}[\vec{\varphi}]=\int D\mu(\vec{\zeta})\int D\mu(\vec{\theta})\ {\rm e}^{-S_{BT}[\vec{\zeta},\vec{\theta},\vec{\varphi}]},$

(38)

where

$\displaystyle S_{BT}[\vec{\zeta},\vec{\theta},\vec{\varphi}]=S_{\rm lead}[\vec{\zeta}]+S_{\rm isl}[\vec{\theta}]+S_{T}[\vec{\zeta},\vec{\theta},\vec{\varphi}],$

(39)

with the lead action

$\displaystyle S_{\rm lead}[\vec{\zeta}(\tau)]=\int_{0}^{\beta}d\tau\sum_{Jk\sigma}\Big{[}\zeta_{Jk\sigma}^{\ast}(\tau)(\partial_{\tau}-\mu+\epsilon_{Jk\sigma})\zeta_{Jk\sigma}(\tau)\Big{]}\equiv\zeta^{\ast}g_{J}^{-1}\zeta\,,$

(40)

and the island action

$\displaystyle S_{\rm isl}[\vec{\theta}(\tau)]=\int_{0}^{\beta}d\tau\sum_{Iq\sigma}\Big{[}\theta_{Iq\sigma}^{\ast}(\tau)(\partial_{\tau}-\mu+\epsilon_{Iq\sigma})\theta_{Iq\sigma}(\tau)\Big{]}\equiv\theta^{\ast}G_{I}^{-1}\theta\,,$

(41)

where $\partial_{\tau}\equiv({\vec{x}_{j}-\vec{x}_{j-1}})/{\Delta_{j}}$ for any variable. We have defined the shorthand notations including the integration over imaginary-time in the vector multiplication. The inverse of the Green’s functions of electrons in the lead $(J)$ and the island $(I)$ are denoted by $g_{J}^{-1}$ and $G_{I}^{-1}$, respectively. We will discuss the Green’s functions in more detail in Section 2.5. The tunneling term reads

	$\displaystyle S_{T}[\vec{\zeta},\vec{\theta},\vec{\varphi}]$	$\displaystyle=\int_{0}^{\beta}d\tau\Bigg{[}\sum_{kq\sigma}\theta_{1k\sigma}^{\ast}t_{1Skq\sigma}\,{\rm e}^{-i\varphi_{1}}\zeta_{Sq\sigma}+\sum_{I=2}^{N}\sum_{kq\sigma}\theta_{Ik\sigma}^{\ast}t_{II-1kq\sigma}\,{\rm e}^{-i(\varphi_{I}-\varphi_{I-1})}\theta_{I-1q\sigma}$
		$\displaystyle+\sum_{kq\sigma}\zeta_{Dk\sigma}^{\ast}t_{DNkq\sigma}\,{\rm e}^{i\varphi_{N}}\theta_{Nq\sigma}+\hbox{H.c.}\Bigg{]}.$		(42)

Inserting the result in Eq.(38) into the representation of the partition function in Eq.(31), we obtain a path integral representation of the partition function which will be the basis for the analysis in the section following. We close this section with a summary of the path integral formulation for the partition function

$$Z[\vec{\zeta},\vec{\theta},\vec{\varphi}]=\mathcal{N}\prod\limits_{I=1}^{N}\bigg{[}\sum_{k_{I}=-\infty}^{\infty}\int\limits_{\varphi_{I}(0)}^{\varphi_{I}(\beta)+2\pi k_{I}}D[\varphi_{I}(\tau)]\bigg{]}\int D\mu(\vec{\zeta})\int D\mu(\vec{\theta})\ {\rm e}^{-S[\vec{\zeta},\vec{\theta},\vec{\varphi}]}\,,$$

(43)

where

$$S[\vec{\zeta},\vec{\theta},\vec{\varphi}]=S_{C}[\vec{\varphi}]+S_{BT}[\vec{\zeta},\vec{\theta},\vec{\varphi}]\,.$$

(44)

Introducing vectors $\vec{\zeta}_{J}=(\ldots,\zeta_{Jk\sigma},\ldots)^{T}$, $\vec{\theta}_{I}=(\ldots,\theta_{Ik\sigma},\ldots)^{T}$ and combining Eqs.(39)–(2.4), we found that

	$\displaystyle S_{BT}[\vec{\zeta},\vec{\theta},\vec{\varphi}]$	$\displaystyle=$	$\displaystyle S_{\rm lead}[\vec{\zeta}]+S_{\rm isl}[\vec{\theta}]+\int_{0}^{\beta}d\tau\Big{[}\theta_{1}^{*}\,\Lambda_{1S}\zeta_{S}$
			$\displaystyle+\sum_{I=2}^{N}(\theta_{I}^{*}\,\Lambda_{II-1}\theta_{I-1})+\zeta_{D}^{*}\,\Lambda_{DN}\theta_{N}+\hbox{H.c.}\Big{]},$

where we have introduced the tunneling matrices

	$\displaystyle{{\Lambda}_{1S}}$	$\displaystyle=$	$\displaystyle{{t}_{1Skq\sigma}}{e}^{-i{{\varphi}_{1}}}$		(46)
	$\displaystyle{{\Lambda}_{II-1}}$	$\displaystyle=$	$\displaystyle{{t}_{II-1kq\sigma}}{{e}^{-i\left({{\varphi}_{I}}-{{\varphi}_{I-1}}\right)}}$
	$\displaystyle{{\Lambda}_{DN}}$	$\displaystyle=$	$\displaystyle{{t}_{DNkq\sigma}}{{e}^{i{{\varphi}_{N}}}},$

and used $\varphi_{S}=\varphi_{D}=0$.

The Green’s functions of the lead electrons are defined by

$\displaystyle(\partial_{\tau}-\mu+\epsilon_{Jk\sigma})g_{Jk\sigma}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime}),$

(47)

and

$\displaystyle g_{Jk\sigma}(0,\tau^{\prime})=-g_{Jk\sigma}(\beta,\tau^{\prime}),$

(48)

where $\mu$ denotes the chemical potential of the lead $J$. Correspondingly, for the island electrons we have

$\displaystyle(\partial_{\tau}-\mu+\epsilon_{Iq\sigma})G_{Iq\sigma}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime})$

(49)

and

$\displaystyle G_{Iq\sigma}(0,\tau^{\prime})=-G_{Iq\sigma}(\beta,\tau^{\prime}),$

(50)

where $\mu$ also denotes the chemical potential of the island $I$. The solution of the inhomogeneous differential equation in Eq.(49) with the boundary condition in Eq.(50) reads [21]

$\displaystyle G_{Iq\sigma}(\tau,\tau^{\prime})=\left\{\begin{array}[]{ccc}{\displaystyle\frac{\exp\left[-\epsilon_{Iq\sigma}(\tau-\tau^{\prime})\right]}{1+\exp\left(-\beta\epsilon_{Iq\sigma}\right)}}&\hbox{for}&\tau>\tau^{\prime}\,,\\ &&\\ -{\displaystyle\frac{\exp\left[-\epsilon_{Iq\sigma}(\tau-\tau^{\prime})\right]}{1+\exp\left(\beta\epsilon_{Iq\sigma}\right)}}&\hbox{for}&\tau<\tau^{\prime}\,,\end{array}\right.$

(54)

with the analogous expression for $g_{Jk\sigma}(\tau,\tau^{\prime})$. For convenience, hereafter, we have measured all energies from the chemical potential of the leads and islands.

The action in Eq.(43) is quadratic in the Grassmann fields $\vec{\zeta}$ and $\vec{\theta}$ as can be seen explicitly form Eq.(2.4). Hence, the corresponding path integral is Gaussian and can be performed analytically. Thus, we can integrate over the quasi-particle reservoirs in three steps as in the following. The first step, using the general Gaussian’s integral formula as [21]

$\displaystyle\int\prod\limits_{i=1}^{N}{\frac{d{{x}_{i}}^{*}\,d{{x}_{i}}}{2\pi i}}\,{{e}^{-x_{i}^{*}\,{{H}_{ij\,}}{{x}_{j}}+\,J_{i}^{*}{{x}_{i}}+\,x_{i}^{*}{{J}_{i}}}}\,\,\,\,={{\left(\det\left(H\right)\right)}^{-1}}{{e}^{\,J_{i}^{*}\,H_{ij}^{-1}\,{{J}_{i}}}},$

(55)

we can integrate over the fermion in the source electrode in Eq.(43) as

	$\displaystyle\int{D\mu\left({{\zeta}_{S}}\right)}\,{{e}^{-S[\vec{\zeta},\vec{\theta},\vec{\varphi}]}}$	$\displaystyle=$	$\displaystyle\int{D\mu\left({{\zeta}_{S}}\right)}\,{{e}^{-\left(\zeta_{S}^{*}g_{S}^{-1}{{\zeta}_{S}}+\theta_{1}^{*}\Lambda_{1S}^{*}{{\zeta}_{S}}+{{\zeta}^{*}}_{S}\Lambda_{1S}\theta_{1}\right)}}$
		$\displaystyle=$	$\displaystyle Z_{S}\,e^{\theta_{1}^{*}\Lambda_{1S}^{*}g_{S}\Lambda_{1S}\theta_{1}},$

where ${{Z}_{S}}=\det\left(g_{S}^{-1}\right)$. The second step, we integrate over ${{\theta}_{1}}$ as

	$\displaystyle\int{D\mu\left({{\theta}_{1}}\right)}\,{{e}^{-S[\vec{\zeta},\vec{\theta},\vec{\varphi}]}}$	$\displaystyle=$	$\displaystyle\int{D\mu\left({{\theta}_{1}}\right)}\,{{e}^{-\left(\theta_{1}^{*}G_{1}^{-1}{{\theta}_{1}}+\theta_{2}^{*}\Lambda_{21}^{*}{{\theta}_{1}}+{{\theta}^{*}}_{1}\Lambda_{21}\theta_{2}\right)+\theta_{1}^{*}\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}\theta_{1}}}$
		$\displaystyle=$	$\displaystyle\det\left(G_{1}^{-1}-\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}\right)\,{{e}^{\theta_{2}^{*}\Lambda_{21}^{*}\tilde{G}_{1}\Lambda_{21}\theta_{2}}}$
		$\displaystyle=$	$\displaystyle{{Z}_{1}}\det\left(1-G_{1}\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}\right)\,{{e}^{\theta_{2}^{*}\Lambda_{21}^{*}\tilde{G}_{1}\Lambda_{21}\theta_{2}}},$

where ${{Z}_{1}}=\det\left(G_{1}^{-1}\right)$ and $\tilde{G}_{1}^{-1}=G_{1}^{-1}-\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}$. We have used the identity matrix properties $G_{1}^{-1}G_{1}=\mathbb{I}$, and the determinant matrix property, respectively. In the same way, one integrates over $\theta_{I}$ variable and then obtains

	$\displaystyle\int{D\mu({\theta}_{I}})\,e^{-S[\varphi,\Lambda]}$	$\displaystyle=$	$\displaystyle Z_{I}\det(1-G_{I}\Lambda_{II-1}^{*}\tilde{G}_{I-1}\Lambda_{II-1})$
			$\displaystyle\times e^{\theta_{I}^{*}\Lambda_{II-1}^{*}\tilde{G}_{I-1}\Lambda_{II-1}\theta_{I-1}},$

where ${{Z}_{I}}=\det\left(G_{I}^{-1}\right)$, and

$\displaystyle\tilde{G}_{I}^{-1}=G_{I}^{-1}-\Lambda_{II-1}^{*}{{\tilde{G}}_{I-1}}\Lambda_{II-1}.$

(59)

The third step, the fermion in drain $\zeta_{D}$ can be integrated as

$\displaystyle\int{D\mu(\zeta_{D})}\,e^{-S[\varphi,\Lambda]}=Z_{D}\det(1-{g_{D}}\Lambda_{DN}^{*}{\tilde{G}}_{N}\Lambda_{DN})\,,$

(60)

where $Z_{D}=\det(g_{D}^{-1})$ and $\tilde{G}_{N}^{-1}=G_{N}^{-1}-\Lambda_{NN-1}^{*}{{\tilde{G}}_{N-1}}\Lambda_{NN-1}$. We may rewrite the partition function in Eq.(38) as

	$\displaystyle Z_{BT}\left[\vec{\varphi}\right]$	$\displaystyle=$	$\displaystyle\mathcal{N}_{BT}\,\det\left(1-G_{1}\Lambda_{1S}^{*}g_{S}\Lambda_{1S}\right)\,$
			$\displaystyle\times\prod\limits_{I=2}^{N}\det\left(1-G_{I}\Lambda_{II-1}^{*}{\tilde{G}_{I-1}}\Lambda_{II-1}\right)$
			$\displaystyle\times\det\left(1-{{g}_{D}}\Lambda_{DN}^{*}{\tilde{G}_{N}}\Lambda_{DN}\right)\,,$

where $\mathcal{N}_{BT}=Z_{S}Z_{D}Z_{1}Z_{2}...Z_{N}$. To simplify future, we rewrite the partition function in Eq.(2.5) in term of the trace of matrix as

	$\displaystyle Z_{BT}\left[\vec{\varphi}\right]$	$\displaystyle=$	$\displaystyle\mathcal{N}_{BT}{{e}^{tr\big{\{}\ln\left(1-G_{1}\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}\right)\big{\}}}}{{e}^{tr\big{\{}\ln\left(1-{{g}_{D}}\Lambda_{DN}^{*}{{\tilde{G}}}_{N}\Lambda_{DN}\right)\big{\}}}}$
			$\displaystyle\times{{e}^{\sum\limits_{I=2}^{N}{tr\big{\{}\ln\left(1-G_{I}\Lambda_{II-1}^{*}{{{\tilde{G}}}_{I-1}}\Lambda_{II-1}\right)\big{\}}}}}$

where we have used the matrix property as

$\displaystyle\det A={{e}^{tr\left\{\ln A\right\}}}.$

(63)

We can therefore rewrite the partition function of the system in term of the effective action as

$\displaystyle Z[\vec{\varphi}]=\mathcal{N}_{sys}\prod\limits_{I=1}^{N}\bigg{[}\sum_{k_{I}=-\infty}^{\infty}\int\limits_{\varphi_{I}(0)}^{\varphi_{I}(\beta)+2\pi k_{I}}D[\varphi_{I}(\tau)]\bigg{]}{{e}^{-{{S}_{\text{eff}}\left[\vec{\varphi},\Lambda\right]}}},$

(64)

where $\mathcal{N}_{sys}=\mathcal{N}\mathcal{N}_{BT}$ and $S_{\text{eff}}\left[\vec{\varphi},\Lambda\right]$ is the effective action of the island system

$\displaystyle S_{\text{eff}}\left[\vec{\varphi},\Lambda\right]={{S}_{C}}\left[\vec{\varphi}\right]+{{S}_{T}}\left[\vec{\varphi},\Lambda\right],$

(65)

with the Coulomb action ${{S}_{C}}\left[\vec{\varphi}\right]$ obtains in Eq.(32) and the tunneling action is expressed as

	$\displaystyle S_{T}\left[\vec{\varphi},\Lambda\right]=$	$\displaystyle-tr\big{\{}\ln\left(1-G_{1}\Lambda_{1S}^{*}{{g}_{S}}\Lambda_{1S}\right)\big{\}}-\sum\limits_{I=2}^{N}{tr\big{\{}\ln\left(1-G_{I}\Lambda_{II-1}^{*}{{{\tilde{G}}}_{I-1}}\Lambda_{II-1}\right)\big{\}}}$
		$\displaystyle-tr\big{\{}\ln\left(1-{{g}_{D}}\Lambda_{DN}^{*}{{{\tilde{G}}}_{N}}\Lambda_{DN}\right)\big{\}}.$		(66)