Magnon-driven quantum dot refrigerators

Magnon-driven quantum dot refrigerators

Physics Letters A 379 (2015) 3054–3058 Contents lists available at ScienceDirect Physics Letters A www.elsevier.com/locate/pla Magnon-driven quantu...

785KB Sizes 0 Downloads 84 Views

Physics Letters A 379 (2015) 3054–3058

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Magnon-driven quantum dot refrigerators Yuan Wang, Chuankun Huang, Tianjun Liao, Jincan Chen ∗ Department of Physics, Xiamen University, Xiamen 361005, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 27 June 2015 Received in revised form 19 October 2015 Accepted 19 October 2015 Available online 23 October 2015 Communicated by A. Eisfeld Keywords: Quantum dot refrigerator Master equation Three-terminal thermoelectric setup Performance evaluation Parametric design

a b s t r a c t A new model of refrigerator consisting of a spin-splitting quantum dot coupled with two ferromagnetic reservoirs and a ferromagnetic insulator is proposed. The rate equation is used to calculate the occupation probabilities of the quantum dot. The expressions of the electron and magnon currents are obtained. The region that the system can work in as a refrigerator is determined. The cooling power and coefficient of performance (COP) of the refrigerator are derived. The influences of the magnetic field, applied voltage, and polarization of two leads on the performance are discussed. The performances of two different magnon-driven quantum dot refrigerators are compared. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The specific low-dimensional structures are attracting much attention due to the advances in techniques for nanostructured materials [1–3]. Of particular interest, since the quantum-dot (QD) refrigerator was presented in 1993 by Edwards et al. [4–7], are the nanothermoelectric setups using the quantum dot system. Thus far, multi-terminal QD setups have been discussed in several references [8–10] because they allow for the crossed flows of the charge current and heat flow. Sánchez and Büttiker proposed a nano-sized structure consisting of two QDs coupled by the Coulomb interaction without particle exchange [8]. EntinWohlman et al. discussed a three-terminal thermoelectric setup composed of a resonant level, two electronic reservoirs, and a phonon source [9]. Li and Jia proposed a particle-exchange heat engine in which three QDs are coupled to two fermionic reservoirs and a bosonic reservoir [10]. The QD coupled to ferromagnetic reservoirs has underlying applications in spintronics, i.e., QD spin valves [11–14]. Strasberg et al. proposed a model of an information driven current through a spin valve with which the Maxwell demon device can be physically realized [15]. Sothmann and Büttiker used two ferromagnetic reservoirs and a ferromagnetic insulator to generate a threeterminal QD setup with only one spin-splitting QD, which can convert part of the heat into an electron current [16]. It gives one type of multi-terminal setup, which can generate a spin-polarized

*

Corresponding author. E-mail address: [email protected] (J. Chen).

http://dx.doi.org/10.1016/j.physleta.2015.10.032 0375-9601/© 2015 Elsevier B.V. All rights reserved.

charge current via a thermal gradient. The electron reservoirs of the conventional thermoelectric devices must be kept at different temperatures and chemical potentials [17–20], whereas the magnon-driven QD setups can be operated between two ferromagnetic metals with identical temperature and used to exploit the heat of a ferromagnetic insulator reservoir [10,16]. Based on the model in Ref. [16], we propose a refrigeration model to cool a bosonic reservoir. Note that when the irreversibility is taken into account, the refrigeration model is not the simple reverse operation of the heat engine model described in Ref. [16]. It includes not only the contribution of some crucial parameters such as the magnetic field, which was not discussed in Ref. [16], but also some applications that may be found in micromechanical systems [21,22] and other fields. In the present paper, the occupation probabilities of the QD are solved by using the model of a magnon-driven quantum dot refrigerator established here and the rate equation. The matter currents are analyzed, and several special cases are discussed. Both the cooling power and coefficient of performance (COP) are optimized by considering the influence of the external magnetic field and applied voltage. The effect of the temperature of the ferromagnetic insulator is considered. The performances of two magnon-driven quantum dot refrigerators with different cooling spaces are analyzed and compared. 2. Quantum dot refrigerator with a ferromagnetic insulator Fig. 1(a) shows a refrigeration system consisting of a spinsplitting QD embedded in two ferromagnetic metallic leads at temperatures T L and T R and chemical potentials μ L and μ R and a

Y. Wang et al. / Physics Letters A 379 (2015) 3054–3058

P0 =

1

Ω 1

P↑ =

Ω

3055

(Γ↑→0 Γ↓→0 + Γ↑→↓ Γ↓→0 + Γ↓→↑ Γ↑→0 ),

(3)

(Γ0→↑ Γ↓→↑ + Γ0→↓ Γ↓→↑ + Γ↓→0 Γ0→↑ ),

(4)

(Γ0→↓ Γ↑→↓ + Γ0→↑ Γ↑→↓ + Γ↑→0 Γ0→↓ ),

(5)

and

1

P↓ =

Fig. 1. (a) The energy diagram of the refrigeration system consisting of a spinsplitting QD embedded into two ferromagnetic metallic leads and an additional ferromagnetic insulator. (b) The transition processes among different states.

ferromagnetic insulator at temperature T B . The two splitting levels are denoted as ε↑ and ε↓ . In Fig. 1(a), γr σ is the spin-dependent tunnel coupling strength between the ferromagnetic metallic lead r (r = L , R ) and the QD with spin σ . The density of state ρr σ of the ferromagnetic metallic lead is spin-dependent. In this model, a polarization p r is introduced, defined as p r = (ρr ↑ − ρr ↓ )/(ρr ↑ + ρr ↓ ). Three special cases, i.e., p r = ±1 and p r = 0, represent completely polarized and unpolarized leads, respectively. The QD in the system is described by the following one-site Hubbard Hamiltonian [16,23,24]



H dot =

σ

εσ d†σ dσ + Ud†↑d↑d†↓ d↓ ,

(1)

where εσ = ε ∓ α /2, the symbol “∓” corresponds to the cases of the energy levels ε↑ and ε↓ , α = μ g g B, B denotes the magnetic †

field, μ g is the Bohr magneton, g is the Lande factor, dσ (dσ ) is the creation (annihilation) operator, and U is the Coulomb energy resulting from the double occupancy of the QD, which is assumed to be infinitely large in the following discussion. The Boltzmann constant k B and electron charge e are set as 1 in the following discussion. The energy exchange between the QD and ferromagnetic metals occurs via electron tunneling, and the occupation of the electron states in the reservoir is given by the Fermi distribution. However, the coupling between the QD and the ferromagnetic insulator described by the Bose distribution is achieved through an exchange interaction. In contrast to the tunnel coupling strengths, the coupling strength γ B between the QD and the insulator will have, in general, nontrivial energy dependence. However, as for the following discussion, this energy dependence is irrelevant [10,16]. The tunnel coupling strength between the ferromagnetic metal1+ p lic lead r and the QD can be expressed as γr ↑ = γr 2 r and 1− p r r 2

γr ↓ = γ

, where γr = γr ↑ + γr ↓ . In sequential tunneling approximation, continuous tunneling of a single electron is defined and the broadening of energy levels can be neglected. Moreover, the collinear magnetizations are taken into account [16]. According to the model illustrated in Fig. 1(b), Eq. (1) results in three quantum states |i  (i = 0, ↑, ↓). The occupation probability of finding the system in a state at time t is denoted as P i (t ) and the energy level is denoted by εi . The evolution of the occupation probability can be written as the rate equation [5,16,17, 25–28]

dP i (t ) dt

=



 Γi  →i P i  (t ) − Γi →i  P i (t ) ,

Ω

where Ω is the normalizing factor. In the symmetric condition, i.e., p L = − p R = p, γ L = γ R = γ B = γ , and μL = −μ R = V /2, the average matter flux entering the QD from the right reservoir can be obtained as

   R↓ R↑ R IM = I M + I M = γ R ↓ P 0 f R− − P ↑ 1 − f R−    + γ R ↑ P 0 f R+ − P ↓ 1 − f R+ ,

where

I rMσ

(6)

is the spin-resolved electron current of lead r with spin

σ and f r∓ = [1 + exp( ε∓α /T2r −μr )]−1 . The magnitude of the mat-

ter flux entering the QD from the right reservoir is equal to that entering the left reservoir from the QD, but their directions are opR L = −I M . Below, we will focus on the case of p > 0. posite, i.e., I M In the absence of a bias voltage, spin-up electrons preferably leave the left reservoir and spin-down ones preferably enter the right one. Although there is no electron exchange between the ferromagnetic insulator and the QD, a magnon current exists due to the exchange interaction between the insulator and the QD. The magnon current can be analytically obtained as B IM =

γ 2

nP ↑ −

γ 2

(1 + n) P ↓ ,

(7)

where n = [exp(α / T B ) − 1]−1 . The heat flows extracted from the three reservoirs are written as

J LQ = (ε↑ − μ L )

1+ p 

2 1− p 

2 1− p 

and B JQ =





γ P 0 f R− − P ↑ 1 − f R−

2 1+ p 

+ (ε↓ − μ R )

 

γ P 0 f L+ − P ↓ 1 − f L+ ,

+ (ε↓ − μ L ) R JQ = (ε↑ − μ R )



γ P 0 f L− − P ↑ 1 − f L−

2



(8)

 

γ P 0 f R+ − P ↓ 1 − f R+ ,

 α  γ nP↑ − (1 + n) P ↓ . 2

(9)

(10)

The external input power is given by R P = (μ R − μ L ) I M .

(11)

In the condition of T B < T E ( T L = T R = T E ), the heat is absorbed from the ferromagnetic insulator and released to the two ferromagnetic leads. The cooling space is the ferromagnetic insulator. The cooling power and COP are given by B Q˙ C = J Q

(12)

and

(2)

i

where Γi →i  is the transition rate from state |i  to |i  . In the limit of the weak contact coupling, only one electron is permitted to exchange during the transition process between the QD and the ferromagnetic metallic leads or the exchange of one magnon changes the spin state of the QD and the ferromagnetic insulator. According to Eq. (2), the occupy probabilities in the steady state can be solved as

ηC =

Q˙ C P

.

(13)

When T B > T E , the heat is absorbed from the two ferromagnetic leads and released to the ferromagnetic insulator. The cooling space is changed to be two ferromagnetic leads. Thus, the cooling power is given by R Q˙ C = J LQ + J Q ,

while the expression of the COP is the same as Eq. (13).

(14)

3056

Y. Wang et al. / Physics Letters A 379 (2015) 3054–3058

Fig. 2. Matter and magnon currents as a function of the applied voltage, where T L = T R = T E , T E = 1.2T B , ε = 0, μ L = −μ R = V /2, α = 2T B , and p = 0.9.

Fig. 3. The schematic diagram of spin-resolved currents, where (a) V = V R , (b) V = V B , and (c) V = V f .

3. Performance characteristic analysis First, the case of T B < T E is considered. As a refrigerator, the heat needs to be absorbed from the cold reservoir and released to the hot reservoir. Because the spin-up electrons preferably tunnel between the left reservoir and the QD and the spin-down ones preferably tunnel between the right one and the QD, the direction of the flow of the matter current is from left to right, and consequently, the relation between the chemical potentials of the two leads should satisfy μ L > μ R . By using Eqs. (6) and (7), the curves of electron and magnon currents varying with the applied voltage can be plotted, as shown in Fig. 2, where two voltages V R and V B R↑ R↓ R↑ L↑ R B denote the cases of I M = 0( I M = − I M ) and I M = 0( I M = − I M ), respectively. The schematic diagrams of spin-resolved currents for the two cases are shown in Fig. 3(a) and (b). Fig. 2 shows that the magnon current monotonously increases with V . When V > V B , the magnon current flows out from the insulator and the electrons flow from the left side to the right side. This shows that the heat is absorbed from the insulator and released to the two ferromagnetic leads. A refrigeration process can be achieved. When V is large enough, the distribution of the two leads f L∓ → 1 and f R∓ → 0, which shows that almost all of the electrons flow from the left side to the right side. In such a case, the magnon current achieves R B the maximum value. When V R < V ≤ V B , I M < 0, and I M ≤ 0, the matter current flows from the left to the right reservoir, accompanied by the magnon current flowing to the ferromagnetic insulator. Electrons tunnel from the hot reservoir to the spin-down level. Afterwards, the electrons flip the spin through the exchange interaction with the insulator, resulting in a magnon current flowing to the insulator. Then, the electrons in the spin-up level tunnel back to the hot reservoir. As a result, the heat is absorbed from the hot reservoir and released to the cold reservoir. When V ≤ V R , R B IM ≥ 0 and I M < 0. The electrons flow from the low to the high potential side. The heat flows from the hot to the cold reservoir. The above analyses show clearly that when V ≤ V B , the system cannot work as a refrigerator. Note that the region between V B and V R decreases as p increases. When p = 1, V B = V R .

Fig. 4. Three dimensional graphs of (a) the cooling power and (b) the COP varying with V and α , where the values of other parameters are the same as those used in Fig. 2.

Based on Eqs. (12) and (13), three-dimensional graphs of the cooling power and COP varying with V and α are shown in Fig. 4. It can be proven from Eqs. (12) and (13) that Q˙ C = 0 and P > 0 when V = V B , and P → ∞ and Q˙ C is finite when V → ∞. Fig. 4 shows that the cooling power is a monotonically increasing function of V , while the COP is not a monotonic function of V . Fig. 4 also shows that Q˙ C = 0 and ηC = 0 when α = 0, and Q˙ C → 0 and ηC → 0 when α → ∞. As a result, the optimal values of α exist at which the cooling power and COP attain their maximum values. Moreover, the boundary between the blue and blank spaces denotes the voltage V B , which increases with α . Using Eqs. (12) and (13), one can draw the curves of the cooling power and COP varying with the applied voltage, as shown in Fig. 5, where α has been optimized. It can be seen that the larger the polarization p is, the larger the cooling power and COP. When L B R p = 1, I M = IM = −I M and the heat and matter currents are proportional to each other. If the QD is initially empty, the spin-up electron tunnels into the QD from the left lead. Afterwards, an exchange interaction between the spin-up electron and the magnon of the ferromagnetic insulator results in a magnon being absorbed from the insulator and flipping the spin of the electron. Then, the spin-down electron tunnels out of the QD to the right lead. During this process, the electron is driven by the interplay of the chemical potential and temperature differences, flowing from the left to the right lead. If the value of V is given, the Carnot COP can be achieved under the condition that the value of α satisfies the reversible transporting condition f R+ = 1 − f L− = n. When the case of T B > T E is considered, the cooling space includes two reservoirs at temperatures T L and T R . The matter current should flow from the right to the left side, and the relation between the chemical potential should satisfy μ L < μ R . It is much more complex to determine the working region of such a refrigerator. By using Eqs. (6) and (7), the curves of spin-resolved electron and magnon currents varying with the applied voltage can be plotted, as shown in Fig. 6. It can be clearly seen that four special operating conditions exist: V = 0, V R , V B , and V f , which will be discussed as follows.

Y. Wang et al. / Physics Letters A 379 (2015) 3054–3058

Fig. 5. (a) The cooling power and (b) the COP as a function of the applied voltage for different polarizations, where α has been optimized and the values of other parameters are the same as those used in Fig. 2.

3057

Fig. 7. Three dimensional graphs of (a) the cooling power and (b) the COP varying with V and α , where the values of other parameters are the same as those used in Fig. 6.

(iii) When V = V B < 0, the magnon current is equal to zero, R↑ R↓ B R i.e., I M = 0, while I M > 0 and I M > I M . According to the particle L↑

Fig. 6. The matter and magnon currents as a function of the applied voltage, where T L = T R = T E , T B = 1.2T E , ε = 0, μ L = −μ R = V /2, α = 2T E , and p = 0.9.

f− r

(i) When V = 0, f r∓ = [1 + exp(

=1−

∓α /2 Tr

)]−1 and it is easy to derive

f r+ . The ratio of the spin-resolved electron currents from R↑ R↓ rσ

the lead r ( I M ) can be determined, i.e., I M / I M = −(1 − p )/(1 + p ). A concise relation between the magnon current and the matter B R current can be found, i.e., I M = −pI M , which is in accord with the result in Ref. [16]. (ii) When V = V R < 0, the matter currents of the system are R L zero, i.e., I M = −I M = 0, but the spin-resolved currents are not R↑

R↓

L↑

L↓

equal to zero. In such a case, I M = − I M > 0 and I M = − I M > 0, as shown in Figs. 6 and 3(a). According to the particle conservaL↑ R↑ B = I M + I M is also larger than zero. tion, the magnon current I M With the help of the matter and magnon currents, heat currents can be analyzed. The heat transformed from the external applied power input is absorbed by the cooling space, and the heat from the insulator is also released to the cooling space. This shows that when V R ≤ V < 0, the heat flows from the hot reservoir to the cold reservoir and the system cannot work as a refrigerator.

R↑

L↓

R↓

conservation, one can easily derive − I M = I M and − I M = I M , as shown in Fig. 3(b). Using Eqs. (6), (11), and (14), one can obtain R Q˙ C = − P = (μ L − μ R ) I M . This means that the total heat transformed from the applied power input is absorbed by the cold reservoir; thus, the system cannot work as a refrigerator when VB ≤ V < VR. R B > 0, I M < 0, and the cooling power (iv) When V = V f < 0, I M is equal to zero. In such a case, electrons release heat to the ferromagnetic insulator. When Q˙ C = 0, the heat transformed from the external applied power input is released to the insulator. The relaB R /(2I M ) between the applied voltage and the magtion V f /α = I M netic field can be obtained from Q˙ C = 0. The schematic diagram of spin-resolved currents is shown in Fig. 3(c). When V f < V < V B , the heat transformed from the external applied power input is divided into two parts, with one part released to the cooling space and the other released to the hot reservoir, such that the system cannot work as a refrigerator. The above analyses show that the system can work as a refrigerator only when V < V f , and consequently, V f is defined as the threshold voltage of a refrigerator. Note that the region between V f and V R decreases as p increases. When p = 1, V f = V B = V R . Based on Eqs. (13) and (14), three-dimensional graphs of the cooling power and COP varying with V and α are shown in Fig. 7. It can be seen that both the cooling power and COP have their respective maximum values. When α → 0, the energy levels and chemical potentials are close to each other, leading to a very large matter current. However, the heat current will approach zero. Moreover, the increase in α results in a decrease in n in the insulator. When α is large enough, the magnon current approaches zero. As a result, the optimal value of α exists. There are two boundaries in Fig. 7. To absorb heat from the two ferromagnetic reservoirs and release heat to the insulator, the applied voltage must satisfy two conditions: V < V f and V > V c , which is a threshold voltage. When V = V c , the heat absorbed from either lead is zero, although

3058

Y. Wang et al. / Physics Letters A 379 (2015) 3054–3058

Fig. 8. The curves of the cooling power varying with the COP, the values of other parameters are the same as those used in Fig. 6.

Fig. 10. The curves of the cooling power varying with the COP, where the black and blue curves denote the cases of T B > T E and T B < T E , respectively.

Fig. 9. Variation of the maximum cooling power with the magnon temperature for different polarizations, where the applied voltage V has been optimized and the values of other parameters are the same as those used in Fig. 6.

cussed in detail by using the spin-resolved currents. The case of T B < T E is discussed in detail. The operating region is determined. The Carnot value of the COP is obtained in the case of p = 1. Moreover, the optimal region of the refrigerator in the case of T B > T E is determined, the maximum values of the cooling power and COP are calculated, and the influence of the ferromagnetic insulator temperature is considered. Moreover, the performances of two refrigerators operated in the cases of T B > T E and T B < T E are compared. It is found that the performances of the refrigerator operating at the condition of T B < T E are much better than those at the condition of T B > T E . This shows that the ferromagnetic insulator, rather than the ferromagnetic reservoir, should be chosen as the cooled space of magnon-driven QD refrigerators. These results may provide some guidance for the optimal design and operation of practical magnon-driven quantum dot refrigerators.

a large electron current is obtained. The values of V c can be obtained through the numerical calculation for p < 1 and V c = −α for p = 1. Using Eqs. (13) and (14), we can plot the curves of the cooling power varying with the COP, as shown in Fig. 8. The maximum values of the cooling power and COP are indicated by Q˙ Cm and ηm C , and the corresponding values of α are defined as (αopt ) Q and (αopt )η . The optimally working region of the refrigerator can be found in the negative slope of the curves in Fig. 8. These regions η Q are indicated by Q˙ C ≤ Q˙ C ≤ Q˙ Cm , ηC ≤ ηC ≤ ηm C , and (αopt )η ≤ α ≤ (αopt ) Q . In Fig. 9, the influence of the temperature of the insulator on the maximum cooling power for a given α is considered. It can be seen that the larger the polarization p is, the larger the span of the temperature. In the absence of a temperature gradient, a matter current driven by the chemical potential and magnon flows from the right to the left lead. The heat absorbed from the two leads is released to the insulator. To compare the performances of the two refrigerators mentioned above, the curves of the cooling power of the two refrigerators varying with the COP are plotted, as shown in Fig. 10. It can be seen that the polarization of the ferromagnetic leads has a huge influence on the performances of the two refrigerators. Fig. 10 shows clearly that in the optimal operation regions, both the cooling power and the COP of the refrigerator operated in the case of T B < T E are obviously larger than those in the case of T B > T E . Thus, one should choose the magnon-driven QD refrigerator operated in the case of T B < T E rather than T B > T E . 4. Conclusions A magnon-driven QD refrigerator consisting of a spin-split QD coupled with two ferromagnetic leads and a ferromagnetic insulator has been established. The working conditions have been dis-

Acknowledgements This work has been supported by the National Natural Science Foundation (No. 11175148), People’s Republic of China. References [1] J.P. Small, K.M. Perez, P. Kim, Phys. Rev. Lett. 91 (2003) 256801. [2] C. Yu, L. Shi, Z. Yao, D. Li, A. Majumdar, Nano Lett. 5 (2005) 1842. [3] R. Scheibner, M. König, D. Reuter, A.D. Wieck, C. Gould, H. Buhmann, L.W. Molenkamp, New J. Phys. 10 (2008) 083016. [4] H.L. Edwards, Q. Niu, A.L. de Lozanne, Appl. Phys. Lett. 63 (1993) 1815. [5] M. Esposito, K. Lindenberg, C. Van den Broeck, Europhys. Lett. 85 (2009) 60010. [6] M. Esposito, K. Lindenberg, C. Van den Broeck, Phys. Rev. Lett. 102 (2009) 130602. [7] M. Esposito, R. Kawai, K. Lindenberg, C. Van den Broeck, Phys. Rev. E 81 (2010) 041106. [8] R. Sánchez, M. Büttiker, Phys. Rev. B 83 (2011) 085428. [9] O. Entin-Wohlman, Y. Imry, A. Aharony, Phys. Rev. B 82 (2010) 115314. [10] P. Li, B. Jia, Phys. Rev. E 83 (2011) 062104. [11] M. Braun, J. König, J. Martinek, Phys. Rev. B 70 (2004) 195345. [12] M.M. Deshmukh, D.C. Ralph, Phys. Rev. Lett. 89 (2002) 266803. [13] S. Braig, P.W. Brouwer, Phys. Rev. B 71 (2005) 195324. [14] B. Sothmann, J. König, A. Kadigrobov, Phys. Rev. B 82 (2010) 205314. [15] P. Strasberg, G. Schaller, T. Brandes, C. Jarzynski, Phys. Rev. E 90 (2014) 062107. [16] B. Sothmann, M. Büttiker, Europhys. Lett. 99 (2012) 27001. [17] H. Wang, G. Wu, Phys. Lett. A 376 (2012) 2209. [18] S. Welack, M. Esposito, U. Harbola, S. Mukamel, Phys. Rev. B 77 (2008) 195315. [19] T.E. Humphrey, R. Newbury, R.P. Taylor, H. Linke, Phys. Rev. Lett. 89 (2002) 116801. [20] H. Wang, G. Wu, Y. Fu, D. Chen, J. Appl. Phys. 111 (2012) 094318. [21] S.H. Ouyang, J.Q. You, F. Nori, Phys. Rev. B 79 (2009) 075304. [22] A. Schliesser, P. Del’ Haye, N. Nooshi, K.J. Vahala, T.J. Kippenberg, Phys. Rev. Lett. 97 (2006) 243905. [23] B. Muralidharan, M. Grifoni, Phys. Rev. B 85 (2012) 155423. [24] L. Ye, D. Hou, R. Wang, D. Cao, X. Zheng, Y. Yan, Phys. Rev. B 90 (2014) 165116. [25] E. Bonet, M.M. Deshmukh, D.C. Ralph, Phys. Rev. B 65 (2002) 045317. [26] D.A. Bagrets, Y.V. Nazarov, Phys. Rev. B 67 (2003) 085316. [27] U. Harbola, M. Esposito, S. Mukamel, Phys. Rev. B 74 (2006) 235309. [28] Y. Zhang, J. He, X. He, Y. Xiao, Phys. Scr. 88 (2013) 035002.