Repositorio Institucional de Documentos

Teor´ıa de scattering de n fotones con m qubits en gu´ıas de onda Advisors: Luis Mart´ın Moreno & David Zueco L´ainez Universidad de Zaragoza - M´aster en F´ısica y Tecnolog´ıas F´ısicas Universit´e de Cergy-Pontoise - Master in Theoretical Physics and Applications Curso 2012/2013 − Eduardo S´anchez Burillo June 24, 2013 1 Abstract Light-matter interaction problem is for sure one of the most important topics in physics. Nowadays it is possible to deal experimentally with few photons, where quantum effects arise. From a theoretical point of view, immediately, numerical approximated methods are necessary. During this work, we study the Matrix Product States (MPS) technique and its application to 1D light-matter interaction problems. Contents 1 Introduction 2 2 MPS 3 2.1 Theoretical fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Expectedvalues................................. 7 2.3 Truncationprocedure.............................. 9 2.4 Timeevolution ................................. 10 2.4.1 Imaginary evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Matrix Product Operators (MPO) . . . . . . . . . . . . . . . . . . . . . . . 12 2.6 WritingstatesasMPS ............................. 13 2.6.1 Productstates.............................. 13 2.6.2 Case with several excitations . . . . . . . . . . . . . . . . . . . . . . 13 3 Physical systems 14 3.1 TightBindingModel.............................. 14 3.2 Decay of an excited qubit coupled to a tight binding model of bosons . . . 16 3.3 Scattering of an incident photon in the RWA regime . . . . . . . . . . . . . 18 3.3.1 Onequbit ................................ 18 3.3.2 Electromagnetic Induced Transparency . . . . . . . . . . . . . . . . 21 3.4 Scattering in the ultra-strong regime . . . . . . . . . . . . . . . . . . . . . 21 3.4.1 Groundstate .............................. 21 3.4.2 Scattering................................ 24 4 Summary, conclusions and future prospects 26 A Details about the truncation 27 1 Introduction Light-matter interaction is a highly important topic in physics. At the quantum level, it is a very rich field, since these systems have a lot of potential (or even real) applications in, for example, quantum cryptography, quantum computing, etc. One of the main problems is to be able to manipulate few photons in the laboratory. Then, in spite of the fact that the field is very interesting, it did not attract too much attention from the theoretical physics community because the theory could not be compared to experimental results. 2 Nevertheless, a lot of effort was put during the last decades from the experimental point of view and these systems are nowadays accessible. For example, interaction of a photon with a single molecule [1] or a superconducting qubit [2, 3] can be implemented in a laboratory, observation of interesting effects like Electromagnetic Induced Transparency (EIT) for a single photon was observed [4], or a problem so exotic as light propagating through a complex network [5] has been studied; even systems mimicking this dynamics, like quantum plasmonics ones [6], have been implemented experimentally, and characteristic quantum properties, like entanglement, have been observed. Consequently the theoretical research has increased exponentially. When studying this kind of problems, usually the analytical models are continuous, but as soon as one wants to study more complicated problems, where analytical tools are not valid, numerical methods are necessary and they require a discretisation of the system, e.g., the space is not continuous but it is just a discrete set of points. Then, the problem becomes a many-body physics one on a lattice. A generic Nsites (that is, a lattice with Nconstituents) state is given by: |Ψi= d1,...,dN X i1,...,iN=1 ci1,...,iN|i1, . . . , iNi,(1) where {|ii}dn i=1 is an orthonormal basis for the Hilbert space Hnof the single site n, dn= dim(Hn) and ci1,...,iN∈C. The Hilbert space of the whole system is H=NN n=1 Hn and its dimension is dim(H) = QN n=1 dn. So, note that we need a number of coefficients which increases exponentially with the number of particles. For simplicity, we can consider dn= 2 ∀nand N= 1000. Then, dim(H) = 21000 >10300. We cannot work with such a large number of coefficients, since it is impossible to store them and, if it were possible, the computations would take unapproachable times. So, it is not possible to handle numerically with the exact expressions, but it is necessary to use approximated methods. In this work we have used the MPS technique to study several light-matter interaction problems in 1D. In the chapter II, we introduce and develop the MPS method. In the chapter III, we present the results. Finally, we present some conclusions and future projects. 2 MPS One of the most celebrated approximated methods in the case of 1D systems with short range interactions is that of Matrix Product States (MPS). As we will see, in the beginning it is shown that any many-body state can be written exactly in such a way that a set of matrices define it. After that, the way in which one takes wisely the approximation will be explained. Lastly, we explain how to deal with this kind of states. 3 2.1 Theoretical fundamentals First of all, the set of coefficients ci1,...,iNcan be stored in a matrix Ci1,(i2,...,iN). On the other hand, every matrix n×madmits a singular value decomposition (SVD) [7], that is, a matrix Acan be written as: A=UΣV†,(2) where Uis an unitary n×nmatrix, Σ is a diagonal n×mmatrix with non negative entries called singular values σi, and Vis unitary m×m. This is applied over C: Ci1,(i2,...,iN)= D2 X n2=1 Ui1,n2σn2V∗ (i2,...,iN),n2,(3) with D2:= min{d1, d2d3. . . dN}. Now, we define (Ai1 1)1,n2:= Ui1,n2and c0 n2,i2,...,iN:= σn2V∗ (i2,...,iN),n2. Then, the state can be written as: |Ψi= d1,...,dN X i1,...,iN=1 D2 X n2=1 (Ai1 1)1,n2c0 n2,i2,...,iN|i1, . . . , iNi(4) We set C0 (n2,i2),(i3,...,iN):= c0 n2,i2,...,iNand we take the SVD of C0: C0 (n2,i2),(i3,...,iN)= D3 X n3=1 U0 (n2,i2),n3σ0 n3V0∗ (i3,...,iN),n3,(5) with D3= min{D2d2, d3d4...dN}. As before, we set (Ai2 2)n2,n3:= U(n2,i2),n3and c00 n3,i3,...,iN:= σ0 n3V0∗ (i3,...,iN),n3. Then: |Ψi= d1,...,dN X i1,...,iN=1 D2,D3 X n2,n3=1 (Ai1 1)1,n2(Ai2 2)n2,n3c00 n3,i3,...,iN|i1, . . . , iNi(6) Repeating the same procedure again and again: |Ψi= d1,...,dN X i1,...,iN=1 D2,...,DN X n2,...,nN=1 (Ai1 1)1,n2. . . (AiN−1 N−1)nN−1,nNc(N−1) nN,iN|i1, . . . , iNi(7) Now, defining (AiN N)nN,1:= c(N−1) nN,iN: |Ψi= d1,...,dN X i1,...,iN=1 D2,...,DN X n2,...,nN=1 (Ai1 1)1,n2(Ai2 2)n2,n3. . . (AiN−1 N−1)nN−1,nN(AiN N)nN,1|i1, . . . , iNi(8) Note that the sums over n2, . . . , nNare in fact products of the matrices Ai1 1,Ai2 2, etc. Then, we rewrite the state: |Ψi= d1,...,dN X i1,...,iN=1 Ai1 1. . . AiN N|i1, . . . , iNi(9) 4 We have shown that every ci1,...,iNcan be written as a product of matrices D1× D2, . . . , DN×DN+1, with D1=DN+1(= 1 in this case). This is called a matrix product state. Because of the way in which we have derived the MPS, D1= 1, but, in general, D1can take other values (e.g., when working with periodic boundary conditions). So, considering D1not necessarily equal to 1: |Ψi= d1,...,dN X i1,...,iN=1 D1,D2,...,DN X n1,n2,...,nN=1 (Ai1 1)n1,n2. . . (AiN N)nN,n1|i1, . . . , iNi,(10) but it is exactly the same as before, that is, a product of matrices, except finally we take the trace over the resultant matrix. Then: |Ψi= d1,...,dN X i1,...,iN=1 Tr(Ai1 1. . . AiN L)|i1, . . . , iNi(11) This is the most general form of a MPS. If a given element of Amis (Aim m)nm,nm+1 , imis called the physical index corresponding to m, since it is related to the local Hilbert space, and the remaining ones are the virtual indices. At the moment, this description is exact for a general many-body state, so it is useless in principle. For simplicity, let di=dand Di=D,∀i. Then, one has dD2Nnumbers to encode dNcoefficients, so, Dwill be an exponentially large number if we want the description to be exact. As it can be supposed, the approximated method will consist in a truncation of the matrix sizes Di. The question is why we can give a reliable description by doing so. When looking at the expression (3), D2is the maximum of {d1, d2d3. . . dN}(in principle, d2d3. . . dN), but some of the singular values could vanish. If χ2is the quantity of non-vanishing σi, one could truncate the value of D2to χ2and the description would be still exact. By repeating the same argument again and again, finally one achieves a state of the following type: |Ψi= d1,...,dN X i1,...,iN=1 χ1,χ2,...,χN X n1,n2,...,nN=1 (Ai1 1)n1,n2. . . (AiN N)nN,n1|i1, . . . , iNi(12) Then, if one has states of this type, it is possible to reduce the values of Dito χi and, hopefully, one will not need anymore an exponential quantity of coefficients. What is more, the singular values could show a huge decay such that it would be possible to truncate Diin such a way that the error committed would remain under control. Anyway, at the moment we do not know if the singular values follow this kind of behaviour. In order to shed light on it, it is not very hard to check that χiis the rank 1of the reduced density matrices 2corresponding to the subsystems consisting of the sites {1, . . . , i}and 1Recall: the rank of a matrix is the number of non vanishing eigenvalues. 2A quantum system can be described by a statistical mixture of pure states, which mathematically is represented by a density matrix [8]. 5 {i+1, . . . , N}3if we consider 1D systems [9]. This is related to one of the most surprising properties of quantum systems: entanglement. A system consisting of two subsystems A and Bis said to be entangled if the state is not a tensor product, that is, if it cannot be written in the following way: |Ψi=|ΦAi|ξBi(13) If a system is entangled, there is not a pure state corresponding to each of the subsystems, but they are described by mixed states, that is, non trivial reduced density matrices. For sure, entanglement is one of the most genuine quantum characteristics, since if one has a classical system consisting of several subsystems, there is a definite state describing each one, instead of having statistical combinations. Obviously, if χis the rank of the density matrices, it is an increasing measurement of the entanglement: if it is one, the state is not entangled since the reduced density matrix corresponds to a pure state, whereas if it increases, the subsystems are further of being pure states, so the system becomes more and more entangled. χhas other interesting properties in such a way that it becomes a good entanglement measurement [9]. Then, coming back to the many-body state written as a MPS (12), we define χ= maxi{χi} as the entanglement measurement. It is reasonable to think that the description of a given slightly entangled state will be faithful if one truncates the values of Di. Actually, χis directly related to the Von Neumann entropy, probably the main entanglement measurement for pure states [10]: S:= −Tr(ρlog ρ),(14) where ρis the reduced density matrix corresponding to one of the subsystems (Sdoes not depend on which subsystem one chooses) and the basis of the logarithm uses to be 2 or e. The relation between χand Sis very close: χsets an upper bound for the Von Neumann entropy [9]. Sturns out to be the capital parameter to justify the MPS method: in [11] they prove that if we want to approximate the ground state of a 1D system with short range interactions by using a MPS, the value of D= maxi{Di}is always bounded by the Von Neumann entropy in the sense that the distance between the most optimal MPS and the true ground state is exponentially small if Dincreases polynomially with the number of particles. It makes perfectly sense; what we are saying is that it is possible to simulate the system via classical computation (because we do not need an exponential quantity of resources) instead of needing quantum computation and it happens when entanglement is small enough, that is, when the system is not so far of being classical in some sense. Thus, that is why MPS method works, even for slightly excited states or time evolution. On one hand, it is widely accepted that physical states live on a tiny subspace of the whole Hilbert space. For instance, when considering time evolution of a given state over a time which scales polynomially with the number of particles, the evolved state samples an exponentially small volume in Hilbert space [12]. On the other hand, ground states of 3If one has a system and considers a partition, the reduced density matrices corresponding to the subsystems have the same rank. 6 local Hamiltonians are slightly entangled, in connection to what we said in the previous paragraph. What is more, they follow the area law [13]. This law arose in the context of black holes theory [14], but it turned out to be a highly general law for quantum field theories and condensed matter physics [13]. It says that if we consider the entanglement entropy Slbetween two parts of the system in such a way that the linear size of the boundary between both parts is l,Slis proportional to the area: Sl∼lD−1, where Dis the spatial dimension, instead of increasing with the volume of the partition, as it does for a general state 4. It confirms the fact that the ground states entanglement is small and it is even much smaller than that of random states. Obviously, the ground state belongs to the little subset of physical states, so it may be expected that the other physical states have small entanglement in such a way that χincreases with some polynomial of the number of sites. MPS is highly related to some older methods, like Numerical Renormalization Group (NRG) [15] and Density Matrix Renormalization Group (DMRG) [16], which were very successful (mainly DMRG) for describing ground states of 1D quantum many-body systems with short range interactions, even though those methods were not very justified in the beginning. Those studies showed that the singular values of ground or even slightly excited states decay exponentially; in other words, these states are slightly entangled, which again is in agreement with the previous analysis. It is important to remark that this method is just valid for 1D systems. On one hand, this is related to the area law, since it means that entanglement is smaller when the spatial dimension decreases. On the other hand, note that the construction of a MPS is adapted to the geometry because of the fact that Dimeasures the entanglement between the piece of the chain at the left of i(including i) and the rest, whereas it does not have this meaning if the system has more dimensions. Anyway, there are similar methods well suited to higher dimensional cases [11], like Projected Entangled Pair States (PEPS) method, but they do not give so good results. In order to illustrate all this stuff, dN numbers are needed to represent a product state and it agrees with D= 1, which is the minimum possible value for D; i.e., if no entanglement (remember that a product state does not have entanglement), Dtakes the smallest possible value. Finally, we have argued that MPS represent well enough low-energy states of some kind of systems. In the following, we will explain how to compute expected values, consider time evolution, find ground states, write some simple states in its MPS form, etc. 2.2 Expected values When one considers different operations with MPS, it is much easier to consider diagrammatic representations. A general coefficient ci1,...,iN, so the whole state, can be 4Critical systems have logarithmic corrections to the area law. In any case, the entropy is still much smaller than that of random states. 7 represented as a network, as it can be looked at the figure (1). There, the boxes are the tensors, the legs pointing up correspond to the physical indices and the others are the virtual ones. The reader can look at other examples of such representations in [17] and [18]. Anyway, it is not so hard to develop the relations by considering the general formula of a MPS, but having the graphical representations is useful to develop a deeper intuition. The coefficients of the state in the chosen basis are: ci1,...,iN= Tr(Ai1 1. . . AiN N) = (Ai1 1)k1,k2(Ai2 2)k2,k3. . . (AiN N)kN,k1,(15) that is, we have Ntensors with three indices and we contract the virtual ones. Note that here the repeated indices are summed; from now on, we will consider this convention. Looking again at the figure (1), one can realise that a contraction corresponds to a link between the virtual indices of neighbours sites. Figure 1: A general coefficient of a state with free boundary conditions is a set of boxes with three legs, except the first and the last ones, where the virtual legs are contracted. This image, and the others of this section, were taken from [18]. Now, if one computes, for example, the square of the norm, it is: hΨ|Ψi=c∗ i1,...,iNci1,...,iN= (Ai1 1)∗ k1,k2(Ai2 2)∗ k2,k3. . . (AiN N)∗ kN,k1(Ai1 1)l1,l2(Ai2 2)l2,l3. . . (AiN N)lN,l1 (16) So, it is just a contraction over the physical indices. The graphical representation is in the figure (2). Note that the physical indices are contracted. Then, we see that it is not necessary to compute all the coefficients, but it is possible to compute things like the norm by processing directly the tensors. Figure 2: When considering the square of the norm, or, more generally, a scalar product, one has to take the graphical representation of both states and contract the physical indices. On the other hand, usually an operator can be written as a sum of products of local operators, that is, a sum of objects like this: 8 O=O1⊗O2⊗···⊗ON,(17) where Onis an operator acting just on the site n. If Oin,jn n:= hin|On|jniand (En)i,j,k,l := (Ain n)∗ i,kOin,jn n(Ajn n)j,l, where nis not summed, the expected value of such a product of local operators is: hOi= (Ai1 1)∗ k1,k2Oi1,j1 1(Aj1 1)l1,l2(Ai2 2)∗ k2,k3Oi2,j2 2(Aj2 2)l2,l3. . . (AiN N)∗ kN,k1OiN,jN N(AjN N)lN,l1(18) = (E1)k1,l1,k2,l2(E2)k2,l2,k3,l3. . . (EN)kN,lN,k1,l1(19) Then, it is possible to compute the expected value of the typical operators without having to obtain explicitly the coefficients, just via the tensors. Now, the graphical representation is in (3). Note that each red box corresponds to one of the Oiand has two physical indices, which are contracted with the ones of |Ψiand hΨ|. Finally, this scheme can not be very efficient. In such a case, as we will explain later, an operator can be written in such a way (its MPO representation) that the expected value can be obtained directly, without having to decompose the operator in a sum of operators like (17). Figure 3: In this case, one contracts the physical indices of the state with those of the local operators. 2.3 Truncation procedure One of the keys of the MPS technique is the truncation, that is, the approximation of a given state by a new one with smaller tensors. Theoretically, the problem is just to minimise the distance between the initial state |Ψi(with tensors An) and a generic one with smaller tensors |Φi(with tensors Bn), that is, we want to minimise: d2(|Ψi,|Φi) = (hΨ|−hΦ|)(|Ψi−|Φi) = 1 + hΦ|Φi−2Re(hΨ|Φi) (20) Then, one has to compute hΦ|Φiand hΨ|Φiand both quantities are calculable. The question is how to minimise such a function in terms of the tensors of |Φi. A possibility is to minimise it iteratively, that is, we start by the first site of the lattice. Then, all the tensors remain constant, except the one corresponding to the first site and the function is minimised with respect to B1; after that, the tensors B1,B3,B4, etc. remain constant and the function is minimised with respect to B2, and so on. When minimising with respect to Bn, if vΨ nand vΦ nare the vectorisations of Anand Bnrespectively, it can be shown (appendix A) that we have to minimise the following quadratic form: f(An, Bn) = 1 + (vΦ n)†E1vΦ n−2Re (vΨ n)†E2vΦ n,(21) 9 0 20 40 60 80 100 1.0 1.5 2.0 2.5 3.0 t <N> Figure 5: Total number of excitations 3.2 Decay of an excited qubit coupled to a tight binding model of bosons After that, we applied the MPS method to a qubit coupled to a tight binding chain of bosons with dipole interaction. So, the Hamiltonian is: H= L X n=1 a† nan−J L−1 X n=1 (a† nan+1 +a† n+1an)+Ωσ+σ−+gσx(a† L0+aL0),(47) where L0= (L+ 1)/2 (Lis taken odd), Ω is the energy splitting of the qubit, gis the coupling constant between the qubit and the guide, σxis the Pauli matrix xand σ±are the ladder operators of the qubit. The interacting term is a point-like dipole interaction: first, the term aL0+a† L0is the electric field in the Coulomb gauge [21]. On the other hand, we suppose that the square modulus of the wave functions corresponding to the states of the qubit ψ0(r) := hr|0iand ψ1(r) := hr|1iare even functions of r, as well as the wave functions can be chosen as real ones. Taking into account that the dipole operator is proportional to r, obviously its diagonal elements will vanish, whereas those non-diagonal will be equal. So, the matrix corresponding to the dipole operator is proportional to σx. The used basis is the same as before (43), except for the site coupled to the qubit, where besides the usual states, we have to add those where the qubit is excited, that is, the local basis is (a† L0)iL0(σ+)n|0i/√iL0!, with n= 0 or 1. We consider the initial condition |Ψ(0)i=σ+|0i, that is, there are not bosons and the qubit is excited. Even though hΨ(0)|N|Ψ(0)i= 1, this dynamics does not preserve the number of excitations, so in principle, it is necessary to consider the whole Hilbert space for the bosons. It would be an unapproachable problem, since the dimension of the space is infinite. However, for the parameters chosen here, it is more than enough with d= 6 (with except d= 12 for L0, since the local Hilbert space is the one corresponding to the boson and the qubit). 16 Now, the local Hamiltonians are: hn=a† nan−J(a† nan+1 +ana† n+1) (n6=L, L0) (48) hL0=a† L0aL0−J(a† L0aL0+1 +aL0a† L0+1)+Ωσ+σ−+gσx(a† L0+aL0) (49) hL=a† LaL(50) We fix the size of the system to L= 21 (so L0= 11). We take = 1, so the hopping factor is fixed to J= 1/π, since the dispersion relation becomes ωk= 2k/π. For the qubit we take Ω = 1 and for the interaction g= 0.05 and g= 0.5. Other details of the simulations are: ∆t= 0.01 and we consider the evolution until t= 10 and the tolerance is 10−6, but the matrices are now allowed to increase, since it is necessary in the case of large gfor the chosen tolerance. The case with g= 0.05 is shown in the figure 6. The gray curve corresponds to the number of excitations and it remains almost equal to 1 (the corrections are at most around 10−4). In such a case, when g/Ω is small enough, the rotating wave approximation (RWA) [22] works well enough; it consists in replacing σx(a+a†) by σ+a+σ−a†(it is said that we are neglecting the counter rotating terms, σ−a+σ+a†; in some sense, in this regime, these induce short time dynamics compared to the rotating terms). The number of excitations is a conserved quantity in this regime, which is in agreement with our results. On the other hand, the black curve corresponds to the population of the excited state of the qubit calculated with MPS. They have been compared to data obtained via exact diagonalisation by using Mathematica and the overlap between both data is total, as it is expected, since the RWA works really well in this case. The case with g= 0.50 is plotted in 7. Again, the gray curve is the number of excitations and it does not remain constant, but it increases. So it is clear that the RWA does not work anymore when g/Ω is large enough. The black one is the population of the qubit computed with MPS. The agreement with [23] seems to be clear and the difference with the RWA results are obvious (not shown). In addition, as well as some of the excitations propagate, others stay close to the qubit, whereas in the case of low g/Ω, all the excitations propagate, nothing remains joined to the qubit (not shown). 17 0 2 4 6 8 10 0.90 0.92 0.94 0.96 0.98 1.00 t N È<eÈpsi>Ȳ Figure 6: Number of excitations (gray) and qubit population (black) for g= 0.05 0 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 t N È<eÈpsi>Ȳ Figure 7: Number of excitations (gray) and qubit population (black) for g= 0.50 3.3 Scattering of an incident photon in the RWA regime Now, we study the scattering problem in the RWA regime through different systems. First, we consider that the scatterer is a qubit and after that we shall deal with the case where there are 2 qubits in order to look for EIT. 3.3.1 One qubit The Hamiltonian for a single qubit coupled to the chain is: H= L X n=1 a† nan−J L−1 X n=1 (a† nan+1 +ana† n+1)+Ωσ+σ−+g(σ−a† L0+σ+aL0) (51) 18 As it was said, the smaller is g/Ω, the better is the RWA. In addition, when considering this kind of problems, the closer is the energy of the incident wave packet to the qubit splitting (Ω), the finer is the approximation. We choose =Ω=1,J= 1/π,g= 0.50 8,L= 170 and L0= 90. Here, the error remains under 10−12, because the system stays in the one particle subspace. We take ∆t= 0.01 and the final time is 230. As initial state we considered the qubit in its ground state and an incident Gaussian wave packet: |Ψ(0)i=X n cna† n|0i, cn=Nexp −(n−n0)2 2σ2+ik0n,(52) where n0= 30, σ= 8 and k0ranges along the set of values {π/2−0.5, π/2− 0.4, . . . , π/2+0.5}, since the interesting physical phenomena appear around the resonance and it happens when the energy corresponding to the wave packet is close to Ω. As ωk= 1 −2/π cos kthe resonance arises at k0=π/2. What we measure is the qubit population as a function of time and the total current flowing between pairs of sites after and before the qubit position. This is just the time integral of the local current [24] jn(t)∝ih(a† n+1an−a† nan+1)i(53) We do not care about the proportionality constant because we are interested in the transmission coefficient, which, for a general incident state |Ψi, is defined as the total transmitted current over the incident one, that is: T(Ψ) := Rdtjn2(t) Rdtjn1(t),(54) where n2is a point beyond the qubit and n1is before the qubit. In the above definition, the integral domain is the time before the photon reflects. On the other hand, it can be computed analytically through the following expression: T(Ψ) = PkkT(k)|Ψ(k)|2 Pkk|Ψ(k)|2,(55) where kruns over π(−1 + 2/L), π(−1 + 4/L),...π,T(k) is the transmission coefficient for a monochromatic wave packet [25] and Ψ(k) := hk|Ψi. The coefficient is shown in the figure (8) and it is compared to the analytical results. As it is seen, both data fit really well, but there are little errors (at most, the relative error achieves 5%) which are due to the time evolution algorithm. The interesting aspect is that the qubit tends to be transparent if the incident energy is far away of Ω, but it reflects almost totally the photon if it is resonant. 8It is true that there are not systems following the Hamiltonian (51); it is just an approximation working at small g/Ω. g= 0.50Ω is clearly out of the RWA regime. However, the intention of this part is to check our method with previous results and the computation time is smaller as gincreases. 19 1.2 1.4 1.6 1.8 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 k THΨL MPS Analytical Figure 8: Transmission factor vs the central momentum for fixed σ= 8. The gray points correspond to those computed by means of MPS technique and the black ones are the analytical results. The qubit population is shown for central momentum π/2−0.5, . . . , π/2 in the figure (9) (for larger values it is not shown because for π/2+0.1 the result almost fits with that for π/2−0.1 and so on). As expected, the qubit excites and after that it relaxes to its ground state. The decay time is bigger if the wave packet is closer to the resonant case, as well as the maximum of the population increases in such a case. 60 80 100 120 140 0.00 0.02 0.04 0.06 0.08 0.10 0.12 t È<eÈpsi>Ȳ k=1.57 k=1.47 k=1.37 k=1.27 k=1.17 k=1.07 Figure 9: Qubit excited state population for several values of the central momentum. As it can be seen, the peak becomes narrower and higher as the wave packet energy is more resonant. 20 3.3.2 Electromagnetic Induced Transparency Now, we consider a photon scattered by 2 different qubits separated by a distance d. In this case, again, if the incident energy is resonant with one of the qubits, then the photon reflects almost totally. However, a new effect appears: if the energy of the photon is in between the energies of both qubits, the pair of qubits becomes transparent. The Hamiltonian is of course the same as before, but now there are 2 qubits instead of 1. The tight binding parameters are the same and those for the qubits are Ω1= 0.6, Ω2= 1.4 and g1=g2= 0.5. The length of the chain is L= 220 and the qubits positions are L1= 109 and L2= 111. The transmission coefficient is shown in the figure (10). The extremal values of the momentum correspond to the resonant cases with both qubits. In the middle, as we knew, a peak close to one appears (it becomes one for the monochromatic case). Again, the results are really close the analytical ones, computed with the transmission factor for this system [25]. 1.0 1.2 1.4 1.6 1.8 2.0 0.0 0.2 0.4 0.6 k THΨL MPS Analytical Figure 10: Transmission factor vs the central momentum for fixed σ= 8. The gray points correspond to those computed by means of MPS technique and the black ones are the analytical results. 3.4 Scattering in the ultra-strong regime Here, we do not neglect the counter-rotating terms, that is, we work beyond the RWA so we take the Hamiltonian (47). First, we search the ground state and after that we study the scattering of a single excitation. 3.4.1 Ground state As it was indicated, when studying scattering experiments, one works at really low temperatures, so the state of the wave guide before putting the excitations is really close to the ground one. Then we have to find it. Note that this is an interesting problem per 21 se; in fact, we have not found any previous work looking for it. We search it by using imaginary time evolution. If gis small enough, the RWA works. In such a case, the ground state is that with no excitations: |0i. However, it is clear that |0iis not an eigenstate of Hif we do not neglect the counter rotating terms. In such a case, the RWA does not work anymore and the system belongs to the ultra-strong coupling regime. Immediately a question arises: are there physical systems working in this regime? The answer is yes. It is known that ordinary systems do not present a so huge coupling, but a lot of effort has been put during the last years and some sophisticated devices, like superconducting qubits, have shown ultra-strong regime [26], besides there are theoretical purposes to achieve even larger coupling constants [27]. For the simulations, we start with small gand the initial state is |0i. After that, we augment gand we use the ground state for the anterior value of gas input state, and so on. We work with L= 20, L0= 10, = 1, J= 1/π, Ω = 1 and g= 0.05,0.10,...,1.00. The algorithm is considered to converge if the difference between subsequent values of the energy is smaller than 10−10 and the difference between subsequent values of the variance hH2i−hHi2is smaller than 10−6. On the other hand, we have again the problem of the infinite Hilbert space on each site, since we are dealing with bosons and there can be any number of excitations, so we have to truncate the local dimension d. We take d= 8 in the qubit position and it decays through the line till it achieves d= 3. A posteriori, we check that this works, since the number of excitations on each site is much smaller than the local value of d. For the matrix sizes, we choose D= 6 in the middle, D= 1 in the limits of the chain and D= 2 in the rest. We show the number of excitations in the figure (11). It achieves its maximum at the qubit position and decays exponentially along the chain. Furthermore, as gincreases, the maximum becomes higher. Finally, as it is seen, it practically vanishes before arriving to the boundaries, so these ground states are the same for larger values of the length (if we took larger values of g, we would have to take larger chains). 22 N_n 0.2 0.4 0.6 0.8 1.0 g 5 10 15 20 n 0.0 0.5 1.0 Figure 11: Population of each cavity as a function of the coupling constant g. Clearly, the peak is an increasing function of g. We show the energy of the ground state vs the coupling constante gin the figure (12). As one can realise, it seems the points lie over a soft curve. If it is true, it implies that there is no crossing between the ground state and the first excited one as gincreases. In [28] they prove that this is true for L= 1 (the Rabi model). Our result would suggest that it is true even for larger values of L. In addition, again in [28], they obtain a behaviour of the energy as a function of gqualitatively like our result. 23 0.0 0.2 0.4 0.6 0.8 1.0 -0.8 -0.6 -0.4 -0.2 0.0 g E_GS Figure 12: Energy of the ground state as a function of g. 3.4.2 Scattering Once we have the ground state, we study scattering. Since we took g= 0.50 in the RWA case, here we choose the same value for the coupling constant in order to compare the results. For the rest of the parameters, we take again the same values as we took in the RWA case (length of the chain, qubit position, width of the wave packet, etc.). Now, for the local dimensions of the Hilbert space, we take the same as we took for the search of the ground state and we impose d= 2 in the rest of the chain, since we expect that, even though the number of excitations is not conserved, the generated excitations remain close to the qubit. On the other hand, we take D= 9 in the middle and it decays until it is D= 2 in the rest of the chain, except for the limits, where D= 1. Because of the limited computational resources, now we need to take ∆t= 0.1. When truncating we stop it if the error is smaller than 10−6or if a whole sweep has been realised. Then, now the results will not be so clean. This is because the error sources are bigger, since the truncation error is higher, the time step is larger and the initial state is not the desired state since there is a little distance between it and the ground state. At the moment we do not have a deep comprehension of the system and we need to gather more data, but we have found interesting behaviours for different values of the incident momentum. We show some of the most significant results for the qubit population in the figure (13). First, note that, before interacting, the qubit population is not zero, since the ground state of the system is not |0i. In addition, note that the curves 24 are now slightly noisy because of the error sources which we mentioned above. Now, we analyse the results. First, for the smallest values of k, the qualitative behaviour is that we observed for the RWA case: the qubit population increases and after that it goes again the its ground state (k= 1.17 in the figure). However, as it increases, the qubit experiences some kind of slow decay (k= 1.52 in the panel). If kis increased, it exhibits again the RWA-like behaviour (k= 1.32. However, after that it again has a new kind of behaviour: it decays slowly but, in addition, it seems that it emits and absorbs again and again part of the excitations (k= 1.72 in the plot). 0 50 100 150 0.08 0.10 0.12 0.14 0.16 0.18 0.20 t È<eÈpsi>Ȳ k=1.72 k=1.52 k=1.32 k=1.17 Figure 13: Qubit population vs time in the ultra-strong coupling regime. For different values of the momentum, qualitatively different behaviours appear. We graph also contour plots with the number of excitations as functions of time and chain position in the figure (14). The left panels (k= 1.17 and k= 1.52) show the normal behaviour: the wave packet splits in transmitted and reflected parts and the rest of the system comes back to the ground state (note that the excitations in the middle are equal before and after the interaction with the wave packet). However, in the top right graphic (k= 1.52) clearly the system does not go immediately to its ground state, but it decays slowly. Finally, in the bottom right panel (k= 1.72) the behaviour is even stranger, since the state in the middle exhibits oscillations, as it happens for the qubit, and in addition it emits subsequent pulses. 25