(*pole representation of the spectral density with prefactor γ and poles in upper half-plane in zlist*)
gam[γ_?NumericQ,zlist_?(VectorQ[#,NumericQ]&),w_?NumericQ]:=γ/(Product[(w-zlist[[i]])(w-Conjugate[zlist[[i]]]),{i,1,Length[zlist]}]);
(*the principal value integral in the mapping*)
intgam[γ_?NumericQ,zlist_?(VectorQ[#,NumericQ]&),w_?NumericQ]:=Re[I (-Sum[γ/(Product[(zlist[[j]]-zlist[[i]])(zlist[[j]]-Conjugate[zlist[[i]]]),{i,1,j-1}]*Product[(zlist[[j]]-zlist[[i]])(zlist[[j]]-Conjugate[zlist[[i]]]),{i,j+1,Length[zlist]}]*(zlist[[j]]-Conjugate[zlist[[j]]])*(zlist[[j]]-w)),{j,1,Length[zlist]}]-1/2gam[γ,zlist,w])];
(*given a coupling strength γ and a vector containing the upper half plane roots and λ and ϵ, this performs the mapping and returns the new γ and the new poles*)
mapping[γ_?NumericQ,zvec_?(VectorQ[#,NumericQ]&),ϵ_?NumericQ,λ_?NumericQ]:=Module[{kvec,rvec},
kvec=Table[I γ/((zvec[[j]]-Conjugate[zvec[[j]]])*Product[(zvec[[j]]-zvec[[i]])(zvec[[j]]-Conjugate[zvec[[i]]]),{i,1,j-1}]Product[(zvec[[j]]-zvec[[i]])(zvec[[j]]-Conjugate[zvec[[i]]]),{i,j+1,Length[zvec]}]),{j,1,Length[zvec]}];
rvec=w/.NSolve[(w-ϵ)Product[(w-zvec[[j]]),{j,1,Length[zvec]}]-Sum[kvec[[j]]*Product[(w-zvec[[i]]),{i,1,j-1}]*Product[(w-zvec[[i]]),{i,j+1,Length[zvec]}],{j,1,Length[zvec]}]==0,w];
rvec=Re[rvec]+I Abs[Im[rvec]]; (*take only upper half plane roots*)
{λ^2 γ,rvec}
]

(*check that the SD is positive*)
Plot[gam[1.,{-2+I,2+I},w],{w,-10,+10},PlotRange->All]

(*just test if it works*)
NIntegrate[gam[1.,{1+I,2+I},w]/(3.-w)*1/(2 Pi),{w,-Infinity,3.,+Infinity},Method->{"PrincipalValue","SingularPointIntegrationRadius"->0.1}]
intgam[1.,{1+I,2+I},3.]

(*the first mapping for a flat SD has not pole decomposition and must be done explicitly, therefore we define a separate distance measre*)
diffmeasure0[Ttg_?NumericQ,ϵtg_?NumericQ,Γ0_?NumericQ,ϵ_?NumericQ,λ_?NumericQ]:=Module[{ret,sptlist},
sptlist=Sort[DeleteDuplicates[{ϵtg-2 Ttg,ϵtg,ϵtg+2 Ttg,ϵ}]];
ret=NIntegrate[(2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*UnitStep[1-((w-ϵtg)/(2 Ttg))^2]-λ^2 Γ0/((w-ϵ)^2+(Γ0/2)^2))^2,
Join[{w,-Infinity},sptlist,{+Infinity}],MaxRecursion->50];
ret
];
(*this defines a distance measure to a semicircle SD characterized by ϵtg and Ttg, allowing variation with respect to ϵ and λ*)
diffmeasure[Ttg_?NumericQ,ϵtg_?NumericQ,γ_?NumericQ,zlist_?(VectorQ[#,NumericQ]&),ϵ_?NumericQ,λ_?NumericQ]:=Module[{ret,sptlist},
sptlist=Sort[DeleteDuplicates[Join[{ϵtg-2 Ttg,ϵtg,ϵtg+2 Ttg},Re[zlist]]]];
ret=NIntegrate[N[(2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*UnitStep[1-((w-ϵtg)/(2 Ttg))^2]-λ^2 gam[γ,zlist,w]/((w-ϵ-intgam[γ,zlist,w])^2+(gam[γ,zlist,w]/2)^2))^2],
Join[{w,-Infinity},sptlist,{+Infinity}],MaxRecursion->50];
ret
];

(*we aim for 2Ttg = 3/4 Γ0 and ϵtg=0*)
(*the optimal choice yields the first Lorentzian spectral density with
(-0.44798874545216505`)^2 Γ0/(w^2+(Γ0/2)^2) *)
NMinimize[diffmeasure0[0.5*0.75,0.0,1.0,0.0,λ],{λ}]
NMinimize[diffmeasure0[0.5*0.75,0.0,1.0,ϵ,λ],{ϵ,λ}]

(*now do the first pole minimization*)
NMinimize[diffmeasure[0.5*0.75,0.0,(-0.44798874545216505`)^2*1.0,{0.5*I},0.0,λ],{λ}]
NMinimize[diffmeasure[0.5*0.75,0.0,(-0.4479887381206968`)^2*1.0,{-1.1052972030642122`*^-10+0.5*I},ϵ,λ],{ϵ,λ}]

(*and the mapping with the optimal parameters*)
mapping[(-0.44798874545216505`)^2*1.0,{0.5*I},0.0,0.3720949697259449`]
mapping[(-0.44798874545216505`)^2*1.0,{0.5*I},-8.043096847226169`*^-9,-0.372094975818078`]

(*yet another minimization, the minimum is already obtained for ϵ->0*)
NMinimize[diffmeasure[0.5*0.75,0.0,0.02778700921459877`,{-0.37174442302717164`+0.25` I,0.3717444230271716` +0.24999999999999994` I},0.0,λ],{λ}]
NMinimize[diffmeasure[0.5*0.75,0.0,0.027787010124485614`,{-0.37174442704872007`+0.24999999729548839` I,0.37174441900562327` +0.25000000270451167` I},ϵ,λ],{ϵ,λ}]

mapping[0.02778700921459877`,{-0.37174442302717164`+0.25` I,0.3717444230271716` +0.24999999999999994` I},0.0,0.3666741914433289`]

NMinimize[diffmeasure[0.5*0.75,0.0,0.003735962351630951`,{-0.5109335299779277`+0.12488220846395447` I,0.` +0.2502355830720911` I,0.5109335299779276` +0.12488220846395434` I},0.0,λ],{λ}]

mapping[0.003735962351630951`,{-0.5109335299779277`+0.12488220846395447` I,0.` +0.2502355830720911` I,0.5109335299779276` +0.12488220846395434` I},0.0,-0.36586118013042696`]

NMinimize[diffmeasure[0.5*0.75,0.0,0.00050007501068037`,{-0.5842712149270853`+0.0681972444081738` I,-0.21196414487137402`+0.1818027555918262` I,0.21196414487137388` +0.1818027555918263` I,0.5842712149270854` +0.06819724440817358` I},0.0,λ],{λ}]

mapping[0.00050007501068037`,{-0.5842712149270853`+0.0681972444081738` I,-0.21196414487137402`+0.1818027555918262` I,0.21196414487137388` +0.1818027555918263` I,0.5842712149270854` +0.06819724440817358` I},0.0,-0.3656755715952492`]

NMinimize[diffmeasure[0.5*0.75,0.0,0.0000668693421556951`,{-0.6272075043681417`+0.04061287557064342` I,-0.34750323387523446`+0.12309139778179518` I,0.` +0.1725914532951229` I,0.34750323387523424` +0.12309139778179515` I,0.6272075043681417` +0.04061287557064332` I},0.0,λ],{λ}]

mapping[0.0000668693421556951`,{-0.6272075043681417`+0.04061287557064342` I,-0.34750323387523446`+0.12309139778179518` I,0.` +0.1725914532951229` I,0.34750323387523424` +0.12309139778179515` I,0.6272075043681417` +0.04061287557064332` I},0.0,-0.365617222323006`]

(*the agreement looks good (this is now for 6 mappings, starting from a constant SD) *)
Plot[{2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*HeavisideTheta[1-((w-ϵtg)/(2 Ttg))^2]/.{ϵtg->0.0,Ttg->0.5*0.75},
gam[0.003735962351630951`,{-0.5109335299779277`+0.12488220846395447` I,0.` +0.2502355830720911` I,0.5109335299779276` +0.12488220846395434` I},w],
gam[8.938823056477509`*^-6,{-0.6541281105210341`+0.025932239839926326` I,-0.43958013516369693`+0.08396410183627712` I,-0.149695618965165`+0.1401036583237967` I,0.14969561896516487` +0.1401036583237968` I,0.43958013516369665` +0.08396410183627676` I,0.6541281105210343` +0.02593223983992638` I},w]},{w,-1,+1}]

(*the first mapping for a flat SD has not pole decomposition and must be done explicitly, therefore we define a separate distance measre*)
diff2measure0[Ttg_?NumericQ,ϵtg_?NumericQ,Γ0_?NumericQ,ϵ_?NumericQ,λ_?NumericQ]:=Module[{ret,sptlist},
sptlist=Sort[DeleteDuplicates[{-ϵtg -2 Ttg,-ϵtg,-ϵtg+2Ttg,ϵtg-2 Ttg,ϵtg,ϵtg+2 Ttg,ϵ}]];
ret=NIntegrate[(2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*UnitStep[1-((w-ϵtg)/(2 Ttg))^2]+2 Ttg Sqrt[1-((w+ϵtg)/(2 Ttg))^2]*UnitStep[1-((w+ϵtg)/(2 Ttg))^2]-λ^2 Γ0/((w-ϵ)^2+(Γ0/2)^2))^2,
Join[{w,-Infinity},sptlist,{+Infinity}],MaxRecursion->50];
ret
];
(*this defines a distance measure to two semicircle SDs characterized by ϵtg and Ttg, allowing variation with respect to ϵ and λ*)
diff2measure[Ttg_?NumericQ,ϵtg_?NumericQ,γ_?NumericQ,zlist_?(VectorQ[#,NumericQ]&),ϵ_?NumericQ,λ_?NumericQ]:=Module[{ret,sptlist},
sptlist=Sort[DeleteDuplicates[Join[{-ϵtg -2 Ttg,-ϵtg,-ϵtg+2Ttg,ϵtg-2 Ttg,ϵtg,ϵtg+2 Ttg},Re[zlist]]]];
ret=NIntegrate[N[(2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*UnitStep[1-((w-ϵtg)/(2 Ttg))^2]+2 Ttg Sqrt[1-((w+ϵtg)/(2 Ttg))^2]*UnitStep[1-((w+ϵtg)/(2 Ttg))^2]-λ^2 gam[γ,zlist,w]/((w-ϵ-intgam[γ,zlist,w])^2+(gam[γ,zlist,w]/2)^2))^2],
Join[{w,-Infinity},sptlist,{+Infinity}],MaxRecursion->50];
ret
];

NMinimize[diff2measure0[0.5*0.75,1.0,1.0,ϵ,λ],{ϵ,λ}]

mapping[1.0,{0.9470760995554037`+I 0.5},0.94707,-0.46840050975031594`]

NMinimize[diff2measure[0.5*0.75,1.0,0.21939903753435583`,{-0.0211727867792757`+0.24999921255078342` I,1.915318886334679` +0.25000078744921655` I},ϵ,λ],{ϵ,λ}]

mapping[0.21939903753435583`,{-0.0211727867792757`+0.24999921255078342` I,1.915318886334679` +0.25000078744921655` I},1.0271298674446205`,0.27370154769073396`]

Plot[{(2 Ttg Sqrt[1-((w-ϵtg)/(2 Ttg))^2]*UnitStep[1-((w-ϵtg)/(2 Ttg))^2]+2 Ttg Sqrt[1-((w+ϵtg)/(2 Ttg))^2]*UnitStep[1-((w+ϵtg)/(2 Ttg))^2])/.{Ttg->0.5*0.75,ϵtg->1.0},gam[0.01643573856275832`,{-0.11448779427902658`+0.20746147911219884` I,1.0138673481171776` +0.09344247347811006` I,2.0218964131618726` +0.199096047409691` I},w]},{w,-2,+2}]

(*Tilmanns parameter:
(Subscript[β, L])=(Subscript[β, R])=β=1
Spektrale Kopplungsstärke Γ=0.02 (wideband)
ϵ=0
Von außen nach innen in der Kette
(Subscript[λ, L])=0.009999959154289148`
(Subscript[λ, 2L])=0.00781394527921633`
(Subscript[λ, 3L])=0.0074508497062710805`
(Subscript[λ, 4L])= 0.007326447541946075`
(Subscript[λ, 5L])=0.007269560988958051`
(Subscript[λ, 6L])=0.007238948047382406`
*)

(0.009999959154289148`/0.02)^2

(*initial SCD: Γ0=1.0, no pole*)
(*first mapping: γ=0.25, pole at +0.5 I*)
mapping[(0.009999959154289148`/0.02)^2,{0.5 I},0,0.00781394527921633`/0.02] 
(*next mappings: as given*)
mapping[0.03816077627390027`,{-0.43301034366239843`+0.25` I,0.4330103436623984` +0.25` I},0,0.0074508497062710805`/0.02]
mapping[0.005296254129782017`,{-0.5672466777480757`+0.1380320898729474` I,0.` +0.22393582025410524` I,0.5672466777480756` +0.13803208987294738` I},0,0.007326447541946075`/0.02]
mapping[0.0007107153788689587`,{-0.6249638132114255`+0.0856770308644717` I,-0.2453426932439562`+0.16432296913552824` I,0.2453426932439562` +0.16432296913552832` I,0.6249638132114254` +0.0856770308644717` I},0,0.007269560988958051`/0.02]
mapping[0.00009389708082947077`,{-0.6548785461957333`+0.057852925704806275` I,-0.3871842149371637`+0.12004400292218843` I,0.` +0.1442061427460106` I,0.38718421493716365` +0.12004400292218852` I,0.6548785461957332` +0.05785292570480626` I},0,0.007238948047382406`/0.02]

(*for sanity check, same for weak coupling...*)
mapping[(0.009999959154289148`)^2*0.02,{0.5 I*0.02},0,0.00781394527921633`] 
mapping[1.2211448407648084`*^-10,{-0.008660206873247967`+0.005` I,0.008660206873247967` +0.005` I},0,0.0074508497062710805]
mapping[6.779205286120979`*^-15,{-0.011344933554961512`+0.0027606417974589473` I,0.` +0.004478716405082105` I,0.011344933554961512` +0.002760641797458947` I},0,0.007326447541946075`]
mapping[3.6388627398090666`*^-19,{-0.012499276264228507`+0.0017135406172894335` I,-0.004906853864879124`+0.0032864593827105653` I,0.004906853864879123` +0.0032864593827105658` I,0.012499276264228509` +0.0017135406172894324` I},0,0.007269560988958051`]
mapping[1.9230122153875605`*^-23,{-0.013097570923914661`+0.0011570585140961273` I,-0.007743684298743275`+0.0024008800584437687` I,0.` +0.0028841228549202126` I,0.007743684298743276` +0.002400880058443767` I,0.013097570923914663` +0.0011570585140961269` I},0,0.007238948047382406`]

Plot[{gam[0.25,{0.5 I},w],
gam[0.03816077627390027`,{-0.43301034366239843`+0.25` I,0.4330103436623984` +0.25` I},w]+1,
gam[0.005296254129782017`,{-0.5672466777480757`+0.1380320898729474` I,0.` +0.22393582025410524` I,0.5672466777480756` +0.13803208987294738` I},w]+2,
gam[0.0007107153788689587`,{-0.6249638132114255`+0.0856770308644717` I,-0.2453426932439562`+0.16432296913552824` I,0.2453426932439562` +0.16432296913552832` I,0.6249638132114254` +0.0856770308644717` I},w]+3,
gam[0.00009389708082947077`,{-0.6548785461957333`+0.057852925704806275` I,-0.3871842149371637`+0.12004400292218843` I,0.` +0.1442061427460106` I,0.38718421493716365` +0.12004400292218852` I,0.6548785461957332` +0.05785292570480626` I},w]+4,
gam[0.000012301073654849794`,{-0.6723305388407864`+0.04150843343787897` I,-0.4747294654205553`+0.08990191865321255` I,-0.17091220073503605`+0.11858964790890857` I,0.17091220073503602` +0.1185896479089086` I,0.4747294654205549` +0.08990191865321231` I,0.6723305388407864` +0.04150843343787902` I},w]+5},{w,-1,+1},PlotRange->All]

Export["test.dat",Table[{w,gam[0.000012301073654849794`,{-0.6723305388407864`+0.04150843343787897` I,-0.4747294654205553`+0.08990191865321255` I,-0.17091220073503605`+0.11858964790890857` I,0.17091220073503602` +0.1185896479089086` I,0.4747294654205549` +0.08990191865321231` I,0.6723305388407864` +0.04150843343787902` I},w]+5},{w,-1,+1,0.001}]]

Plot[gam[0.000012301073654849794`,{-0.6723305388407864`+0.04150843343787897` I,-0.4747294654205553`+0.08990191865321255` I,-0.17091220073503605`+0.11858964790890857` I,0.17091220073503602` +0.1185896479089086` I,0.4747294654205549` +0.08990191865321231` I,0.6723305388407864` +0.04150843343787902` I},w],{w,-1,+1}]
Plot[gam[1.0077039538052946`*^-27,{-0.013446610776815728`+0.000830168668757581` I,-0.009494589308411098`+0.0017980383730642519` I,-0.0034182440147007198`+0.0023717929581781707` I,0.003418244014700719` +0.0023717929581781702` I,0.009494589308411109` +0.001798038373064254` I,0.01344661077681573` +0.0008301686687575801` I},w],{w,-0.02,+0.02}]

trans[γ_,zvec_,w_]:=Re[gam[γ,zvec,w]^2/((w-1/Pi NIntegrate[gam[γ,zvec,w1]/(w-w1),{w1,-Infinity,w,+Infinity},Method->{"PrincipalValue","SingularPointIntegrationRadius"->0.01}])^2+gam[γ,zvec,w]^2)];

ListPlot[Monitor[Table[{w,trans[0.000012301073654849794`,{-0.6723305388407864`+0.04150843343787897` I,-0.4747294654205553`+0.08990191865321255` I,-0.17091220073503605`+0.11858964790890857` I,0.17091220073503602` +0.1185896479089086` I,0.4747294654205549` +0.08990191865321231` I,0.6723305388407864` +0.04150843343787902` I},w]},{w,-1,+1,0.002}],w]]

Export["test.dat",Monitor[Table[{w,trans[0.000012301073654849794`,{-0.6723305388407864`+0.04150843343787897` I,-0.4747294654205553`+0.08990191865321255` I,-0.17091220073503605`+0.11858964790890857` I,0.17091220073503602` +0.1185896479089086` I,0.4747294654205549` +0.08990191865321231` I,0.6723305388407864` +0.04150843343787902` I},w]},{w,-1,+1,0.002}],w]]

ListPlot[Monitor[Table[{w,trans[1.0077039538052946`*^-27,{-0.013446610776815728`+0.000830168668757581` I,-0.009494589308411098`+0.0017980383730642519` I,-0.0034182440147007198`+0.0023717929581781707` I,0.003418244014700719` +0.0023717929581781702` I,0.009494589308411109` +0.001798038373064254` I,0.01344661077681573` +0.0008301686687575801` I},w]},{w,-0.02,+0.02,0.002*0.02}],w]]