Reduced Density-Matrix Approach to Strong Matter-Photon Interaction

We present a first-principles approach to electronic many-body systems strongly coupled to cavity modes in terms of matter–photon one-body reduced density matrices. The theory is fundamentally nonperturbative and thus captures not only the effects of correlated electronic systems but accounts also for strong interactions between matter and photon degrees of freedom. We do so by introducing a higher-dimensional auxiliary system that maps the coupled fermion-boson system to a dressed fermionic problem. This reformulation allows us to overcome many fundamental challenges of density-matrix theory in the context of coupled fermion-boson systems and we can employ conventional reduced density-matrix functional theory developed for purely fermionic systems. We provide results for one-dimensional model systems in real space and show that simple density-matrix approximations are accurate from the weak to the deep-strong coupling regime. This justifies the application of our method to systems that are too complex for exact calculations and we present first results, which show that the influence of the photon field depends sensitively on the details of the electronic structure.


Short Guide to this supporting information
In this supporting information, we present in detail how we performed the convergence study of the numerical examples shown in the paper and especially how we validate our implementation in Sec. 1 to 3. We provide a "HowTo" for such calculations in Sec. 4. In the last section, we examine the issue of the bosonic symmetry in the photon wave function, which is relevant for the argumentation in the main paper.
To perform the calculations, we integrated our routine in the Octopus code. 1 To find the accuracy threshold of our implementation we compared the Octopus results to calculations of a private code "Dynamics" that was used and validated for Ref. S1 by S.E.B. Nielsen. 2 Both codes approximate the polaritonic orbitals on discretized real-space boxes in the x and q coordinate, but they use different boundary conditions, a different level of finite differences for the approximation of the differential operators (fourth order in Octopus vs sixth order in Dynamics), a different method of the grid point evaluation (on-point in Octopus vs. mid-point in Dynamics), and also a different orbital optimization technique. Specifically, we present the convergence studies for the example of the one-dimensional Helium atom inside a cavity (He) that we used among others in the main part of this paper. The He atom is described by the nuclear potential v He (x) = − 2 √ x 2 +1 (we use atomic units throughout) and coupled to one cavity mode with frequency ω and coupling parameter λ. The modified local potential due to the dressed auxiliary construction reads v He (x, q) = v He (x)+ 1 2 (λx) 2 + ω 2 2 q 2 − ω v (x, q), and the modified interaction kernel w (xq, x q ). Additionally, the many-body wave function that is constructed from dressed orbitals exhibits an additional exchange-symmetry between the q-coordinates. In this work, we do not take this symmetry into account, which is a reasonable approximation for the presented calculations. For further details see Sec.
The convergence study and validation of our implementation in Octopus requires several steps, which correspond to the sections of this supporting information. We start with converging the exact dressed many-body ground-state and compare it to the results of Dynamics. Although we are limited to a very small Hilbert space for such exact calculations, we only have to solve a linear eigenvalue problem, so it is easy to converge the results to a very high precision within the accessible parameter range. Thus, we obtain an upper bound for the convergence accuracy from these exact calculations. In the next section, we present the convergence of dressed HF (that we abbreviate by dHF in this supporting information) in all its details and compare the converged ground-state to Dynamics. This allows us to determine the accuracy for dHF calculations, that is (because of its non-linear nature) harder to converge than the exact routine. As a side effect, the comparison of the two codes is a good validation of the correct implementation of the dressed modifications. Since these modifications are exactly the same for dressed RDMFT (that we abbreviate by dRDMFT in this supporting information) and dHF, validating the dHF implementation validates also the corresponding changes for dRDMFT. In the last section, we discuss the convergence of dRDMFT, which requires two different minimization procedures that are interdependent, and thus again is harder to converge than dHF.
All the above sections are organized similarly. We start with the separated problem of the electronic system outside the cavity and converge first the electronic counterparts of the routines, i.e. the electronic exact solver, HF or RDMFT. The purely photonic problem of the separated problem remains the same at all the discussed levels of theory. We therefore discuss it only once in the first section. Then, we analyze the convergence of the dressed S3 theories in the no-coupling (λ = 0) limit, which theoretically means that the electronic and photonic problems are perfectly separated, but solved in only one large calculation. By comparison of only the electronic (photonic) part of the dressed solutions to the results of the purely electronic (photonic) theories, we can measure if there is a decrease in accuracy due to the simultaneous description of electronic and photonic coordinates. We finish the sections with a convergence study for λ > 0.
1 Validation of the exact dressed many-body groundstate We start with the validation of the exact many-body ground-state. In both codes, this is calculated by minimizing directly the energy expression of the full many-body wave function (denoted by Ψ in the main text), discretized on the grid. For the He test system, Ψ = Ψ (x 1 q 1 , x 2 , q 2 ) is four-dimensional. The actual minimization is performed in Octopus by conjugate gradients method, whereas Dynamics makes use of a Lanczos algorithm.
At the beginning of every calculation, we need to find the proper grid, which is defined by the box sizes L x , L q and the spacings dx, dq in both dimensions. For the minimization in Octopus, we used two different convergence criteria, E = 10 −9 and ρ = 10 −8 . The former tests the energy deviations and the latter the integrated absolute value of the density deviations between subsequent iteration steps. 3 These are the criteria already available in Octopus, so we here show that these are sufficient to produce reliable results. Note that we choose E = 1 10 ρ , because ρ is a much stricter criterion. For the box size and spacing convergence, we perform series C = {C 1 , C 2 , . . . }, where C can be L x , L q , dx, or dq in the following. We perform two types of convergence tests for every C. In the first one, we investigate the deviations in the energy between subsequent elements ∆E We denote the corresponding thresholds with E C . The second type of convergence series considers the maximal deviations in the absolute-value of the electronic/photonic part of the polaritonic density ∆ρ for the series L q , dq. We denote the corresponding thresholds with ρ C .
We start with the electronic part of the example (corresponding to the He atom outside the cavity), choose dx = 0.14 and vary L x = {8, 10, 12, . . . }. The ground-state energy E Lx drops with increasing L x (because boundary effects become less important) and we find ∆E Lx < 10 −8 and ∆ρ Lx < 10 −7 for L x ≥ 20. In the following, we choose L x = 20 because the thresholds of ρ Lx = 10 −7 and E Lx = 10 −8 are already stricter than the maximal accuracy of the later non-linear (dHF,dRDMFT) calculations. Next, we perform dx = {0.2, 0.19, 0.18, . . . , 0.06}. We find that ∆E dx also drops with decreasing dx until ∆E dx < 10 −8 for dx = 0.15 and does not decrease any more (as the convergence criteria of the minimization are no more precise.) ∆ρ dx instead decreases (slowly) until the lowest tested value of dx = 0.06. We find ∆ρ dx < 10 −7 already for dx < 0.14, but to reach ∆ρ dx < 10 −8 , we need to decrease the spacing until dx = 0.07. Such a small spacing is numerically unfeasible for larger boxes and thus we will not try to go beyond ρ dx = 10 −7 .
We repeat the same series for the harmonic oscillator system of the q-coordinate with frequency ω = ω res ≈ 0.5535 (resonance with the transition between ground and first excited state of the He atom outside the cavity) and find that the box size is converged with ∆E Lq < 10 −8 and ∆ρ Lq < 10 −8 for L q ≥ 14. The spacing is converged in the energy ∆E dq < 10 −8 and in the density ∆ρ dq < 10 −7 for dq ≤ 0.20.
From these preliminary calculations, we can infer the parameters for the dressed calculations: L x = 20, L q = 14, dx = 0.14, dq = 0.20. To be sure that these parameters are still sufficient for non-zero coupling (λ > 0), we perform another box length series in the four-dimensional space of the exact ground-state, e.g. for λ = 0.1 and ω = ω res . To have a uniform grid distribution, we also set dx = dq, although our preliminary calculations suggest that we could choose a larger dq. Unfortunately, we cannot explore this space completely, but we are limited with L x , L q ≤ 18. 4 However, we can confirm that the energy is converged in the q-direction for L q = 14 and in the x-direction, we have ∆E Lx ≈ 10 −7 for L x = 18.
Finally, we compare the Octopus results with the results from the Dynamics code. For that, we consider the differences in the total energy E OD = |E Dynamics − E Octopus | and the maximum deviations in the polaritonic density ρ OD ≡ max x,q |ρ Octopus (x, q) − ρ Dynamics (x, q)| between the two codes. Due to the optimization of Dynamics, we are limited to a box length of L x = L q = 14 to have a spacing of dx = dq = 0.14. With these parameters, all the energy and density errors in Octopus increase to ∆E Lx,Lq,dx,dq = ∆ρ Lx,Lq,dx,dq ≈ 10 −5 . When we compare the ground-state energies of both codes for this maximum possible mesh, we find E OD ≈ 10 −5 and ρ OD ≈ 10 −5 . So both codes agree on the level of accuracy that we estimated from the Octopus calculations before, although they are quite different from a numerical perspective. We conclude from these calculations that the implementation of the exact many-body routine for dressed two-electron systems 5 in Octopus is reliable and we can use it as benchmark for the dHF and dRDMFT approximations within Octopus.

Validation of the dHF routine
In this section, we present the validation of dHF, which requires several steps. We start with the electronic HF routine and converge the He atom model (outside the cavity) in box size and spacing, where we proceed like in the previous section. Then, we converge a second HF routine (HF basis ) that is implemented in Octopus, which makes use of a basis set. The difference of such a basis-set implementation is that the routine calculates all the integral kernels of the total energy as matrix elements of a chosen basis and then searches for the energy minimum by only varying the corresponding coefficients. This routine can be considerably faster than the standard one of Octopus, because the calculation of the exchange-term, which is numerically much more expensive than all the other parts of the Hamiltonian, needs only to be performed once in the beginning for the matrix elements.
The conjugate gradient algorithm of the standard HF routine of Octopus instead needs to evaluate the Hamiltonian including the exchange term at every iteration step. Such a basisset implementation requires a convergence study with respect to the basis size, which we will explain in detail in the second part of this section together with a comparison between HF 5 Larger systems are in principal also possible, but numerically infeasible in real space.

S7
and HF basis . We conclude the section with the discussion of dHF, which also uses a basis set.
However, the basis-set convergence for dHF is more involved than for the purely electronic HF basis and we explain this in detail. Afterwards, we compare dHF in the no-coupling limit to HF/HF basis and we discuss the comparison of Octopus and Dynamics on the level of dHF, where we follow the same strategy like in the previous section.

HF minimization on the grid
We start with the He atom outside the cavity and calculate the HF ground-state with the standard routine of Octopus, which uses a conjugate gradients algorithm for the orbital optimization. The convergence criteria for the conjugate gradient algorithm are defined exactly as before for the exact many-body ground-state. We set them to E = 10 −9 and ρ = 10 −8 and determine L x and dx. We start with a spacing of dx = 0.1 and find that ∆E Lx ≈ 10 −8 and ∆ρ Lx < 10 −7 for all L x ≥ 20. We choose L x = 20 and perform the spacing series, which is converged in energy with ∆E dx < 10 −8 for dx ≤ 0.15 and in the density with ∆ρ dx < 10 −7 for dx ≤ 0.14. These calculations suggest that all dx ≤ 0.14 are sufficient and we choose dx = 0.1, which is numerically feasible for the test-system, we consider.

HF minimization with a basis set
In this section, we present the convergence of the He test system in the HF basis routine with respect to the basis set. As Octopus is a real-space code, there are no standard quantumchemistry basis sets implemented. This means that if we wish to express the Hamiltonian and the wave function in a basis set, we need to generate one. For all the calculations of this paper, we do this in the same way by performing a preliminary calculation with the independent particle (IP) routine, which solves a simple 1-body Schrödinger equation for every orbital, neglecting interaction. Among the numerically efficient options that Octopus offers, the basis from IP converged fastest with respect to the size of the basis set M . 6 To 6 Note that we used the symbol M in the main text.

S8
allow for enough variational freedom, we calculate besides the occupied orbitals 7 (that form the HF ground state and are denoted by GS) also excited orbitals, which are called extra states (ES) in Octopus. The so generated basis set has the size M = GS + ES.
The routine performs the minimization by representing the coefficients of the basis-set expansion in a matrix and diagonalizing it repeatedly until self-consistence. The corresponding energy convergence criterion E remains the same in HF basis like in HF. But the convergence criterion Λ of the orbital minimization needs to be adapted and there are several possible options. In Octopus, Λ tests the hermiticity of the Lagrange multiplier matrix Λ S2 that guarantees the orthonormality of the natural orbitals. This is not an obvious choice, but it is more general than typical criteria, and thus allows for using it also for RDMFT minimizations in a basis set. In fact, HF basis and RDMFT in Octopus use exactly the same algorithm for the orbital optimization. This method is explained in detail in Ref. S3. As Λ is a considerably stricter criterion than E , it is set as default to Λ = 10 3 · E .
For the convergence of the size of the basis set M or equivalently the parameter ES, We use L x and dx that we determined before in HF and find that ∆E ES,100 and ∆ρ ES,100 decrease with increasing ES and for ES ≥ 40, ∆E ES,100 < 10 −8 and ∆ρ ES,100 ≈ 10 −5 .
Comparing HF with HF basis with these parameters, we find that |E HF − E HF basis | ≈ 10 −8 and max x |ρ HF (x) − ρ HF basis (x)| ≈ 10 −5 . So both methods are consistent.

Validation of dHF
After having properly understood the convergence of the electronic HF methods, we can now turn to dHF. All the calculations shown in the main part of this paper are done with the basis-set type of implementation, which is more involved than in the purely electronic 7 Note that in our case, we have always GS = N 2 , where N is the number of electrons, because we only look at closed-shell systems, which distribute two electrons to every spatial orbital.

S9
HF basis case. Thus, we start in Sec. 2.3.1 with discussing the new issues that enter when one needs to converge a system in the dressed space with respect to the basis set. In Sec.
2.3.2, we present the convergence of dHF with zero-coupling (λ = 0), which means that the electronic and photonic part of the system completely decouple such that we can compare the results of the electronic part to electronic HF. In the last section, 2.3.3, we compare the dHF results from Octopus with Dynamics.

Basis-set convergence in the dressed auxiliary space
In the dressed auxiliary space that we explore with dHF and dRDMFT, the basis-set convergence is not as straightforward as in the HF basis case. We illustrate this additional difficulty for dHF in the no-coupling (λ = 0) case. λ = 0 makes the discussion much simpler, because we have two perfectly separated problems, the atomic system and harmonic oscillators that are just calculated at the same time. However, we use a combined basis that consequently needs to include the appropriate degrees of freedom for both systems. But the composition of the basis set is strongly influenced by the parameters of the system, especially ω, because it is generated by a preliminary calculation. This can be illustrated as follows: For λ = 0, the electronic and photonic part of the system separate. Thus, the coupled (dHF) Hamiltonian is a direct sum of the electronic (photonic) HamiltonianĤ e (Ĥ p ) and the two-dimensional dressed orbitals ψ iα (x, q) can be exactly decomposed in their one- Here, φ i (χ α ) is an eigenfunction with eigenvalue e e i (e p α ) of the electronic (photonic) Hamil-tonianĤ e (Ĥ p ). Consequently, we know that we can calculate the eigenvalues of the dressed orbitals as sum of the uncoupled ones, i.e. e ep iα = e e i + e p α . The basis set for a dHF calculation is then constructed using the ground and the first ES excited orbitals of a preliminary IP calculation. These orbitals are ordered by their eigenvalue e ep iα and for that, the relation between the individual energies of the electron and photon space is crucial. Table S1 shows the decomposition of such a basis with ES = 6 (and thus M = 7.)

S10
We see that slightly more electronic states contribute to the basis, but already 3 out of 7 states describe excited photon contributions. So if we wanted to use that basis for an dHF calculation with for example λ = 0, that trivially needs only ground-state contributions for the photonic coordinate, these 3 states would be entirely unnecessary and waste computational resources. Before we discuss this further, we want to explain how ω influences this distribution between electronic and photonic contributions: If we chose for example a larger ω ≈ 1.3, the first 4 states of the combined basis would only vary in the electronic contribution. If we instead lowered ω, the opposite would happen and the first states would vary in the photonic contribution. A similar kind of argument could be done for all other ingredients of the Hamiltonian of course, but for ω this influence is most directly visible and comparatively strong. A detailed investigation of this issue is beyond the scope of this paper.
Pragmatically, we sometimes adapt ω such that the produced basis is reasonable (see the Beryllium test system in Sec. Numerical Results of the main part of the paper.) At this point, one might ask the question why at all we use the basis-set implementation and the answer remains the same as for the purely electronic case: Although we need large basis sets for the properly converged dHF calculations which makes them numerically very expensive, these computations are still relatively inexpensive compared to calculations with a conjugate gradients algorithm that calculates all integrals on the grid. In the case of the He atom, the dHF calculation with conjugate gradients and ES = 15 takes still 4 times longer than the same calculation with a basis set and ES = 50.
Finally, we want to mention that the statements of this subsection carry over straightforwardly to the coupled (λ > 0) case, we just do not have the direct connection to the decoupled spaces and thus cannot visualize this case so clearly. It is clear that we need more variational freedom for the photonic subspace than in the λ = 0 case but the basis will still depend on the system parameters. One can expect that also for λ > 0 many basis states will be unnecessary. of the combined HamiltonianĤ ep and the decomposition of the combined index: The first eigenenergy e ep 1 = e ep 11 = e e 1 + e p 1 is the sum of the two first uncoupled energies. The second energy e ep 2 = e ep 12 instead is formed from the electronic ground but photonic first excited state, etc. In this example, we see that both contributions are similar, although there are slightly more electronic orbitals included. This tendency continues also for increasing ES. Now we can address the ES-convergence of dHF, which we will present for λ = 0, such that we can compare the electronic part of the converged result to HF afterwards. As we set λ = 0, we know that for the photon component, the IP ground-state orbital χ 1 (q) is already the exact solution, because the photons do not interact directly with each other and HF (or any other level of approximation) will not change that. The dHF ground-state orbital thus is ψ(x, q) = φ HF (x) ⊗ χ 1 (q) with the electronic HF ground-state orbital φ HF (x). As we saw in the last section, when we converged the He test system using the HF basis routine, φ HF (x) can be approximated by an IP basis with ES = 40. So we know that in the no-coupling limit of dHF, we do not need a larger basis set, because Kronecker-multiplying the basis-states of the HF basis -calculation with χ 1 (q) would be exactly sufficient. However, here and in the following, we do not want to choose by hand such a well-adapted state space (which for very strong coupling strength is non-trivial), but instead rely on the polaritonic IP-calculations.
From our considerations before, we expect to need more basis states to reach the same level S12 of accuracy than using HF basis and the exact number will depend on ω. Indeed, for ω = ω res , we need ES = 80 to have ∆E ES,100 < 10 −7 and ∆ρ ES,100 ≈ 10 −5 , where we calculate the electron density of the dHF-calculation by ρ(x) = dq ρ(x, q). When we instead choose a very high value of ω = 5.0, we reach ∆E ES,100 < 10 −7 and ∆ρ ≈ 10 −5 already for ES ≥ 40.
For both calculations, we do not reach ∆E ES,100 < 10 −8 even for higher ES. This is due to the many "unnecessary" states that are taken into account in the minimization. These then only introduce numerical noise without providing useful variational information. We find this confirmed by analyzing the "photonic density" ρ(q) = dx ρ(x, q), which deviates from the correct density ρ ex (q) = 2|χ 1 | 2 (q) although χ 1 (q) is explicitly part of the basis. Adding basis states consequently cannot improve the photonic orbital. For instance, for ω = ω res , the density deviations remain at approximately ∆ρ ES,ES ref ≈ 10 −4 for all ES and ES ref .
Despite being less accurate, we used ω = ω res for the results section in the main part of the paper (Sec. Numerical Results,) because ∆ρ ES ≈ 10 −4 is still two orders of magnitude smaller than the deviations from the exact or the dRDMFT solution. 8 Finally, we want to compare the electronic part of dHF to HF. For the energy comparison, we need to subtract the photon part E p = 2 ω 2 from the total dHF energy E dHF , E dHF,e = E dHF − E p . However, this analytical expression is also not exact, because of the just mentioned error in the photonic part of the orbitals. Consequently, we have deviations for ω = ω res of |E HF − E ω=ωres dHF,e | ≈ 3 · 10 −6 and max x |ρ HF (x) − ρ ω=ωres dHF,e (x)| ≈ 10 −4 , but for |E HF − E ω=5.0 dHF,e | ≈ 4 · 10 −7 and max x |ρ HF (x) − ρ ω=5.0 dHF,e (x)| ≈ 10 −5 .

Validation of dHF with Dynamics
We conclude this section by comparing the results of our implementation of dHF in Octopus with the Dynamics code, which uses an imaginary time propagation algorithm to calculate the HF ground-state. This comparison allows us to validate also the case of λ > 0. We choose ω = ω res and λ = 0.1 for the comparison and perform another ES-convergence, confirming that we need ES = 90 for the same level of accuracy like before in the no-coupling case. For the sake of completeness, we want to mention that for ω = 5.0 we need again a significantly smaller basis set for even better converged results, exactly as we found before in the λ = 0 case. However, already for ω = ω res the level of convergence is an order of magnitude better than the expected deviations between the codes that we estimated in Sec. 1.
So we compare the ground-state for these parameters to the result of Dynamics and find for the energy the expected deviations of E OD ≈ 10 −5 . For the density, we find instead deviations of ρ OD ≈ 10 −3 . This discrepancy is probably due to the different convergence criteria of the two codes. Dynamics tests the eigenvalue equation of the one-body Hamiltonian for a certain subset of all the grid points, which is a much stronger criterion than the one of Octopus, explained in Sec. 1. The influence of these different criteria on a self-consistent calculation are naturally stronger than on the calculation of a linear eigenvalue-problem like the many-body calculation. Still, density errors of the order of 10 −3 are small enough for our purposes. We find similar errors also for other values of λ and conclude that both codes are sufficiently consistent.

Validation of dRDMFT
In this last section we turn to dRDMFT. This is the first implementation at all of this theory, so we cannot validate it with a reference code any more. However, the difference between HF and RDMFT (using the Müller functional) on the implementation-level essentially is in the treatment of the occupation numbers, which are fixed to 2 and 0 in the former case but are allowed to be non-integer for the latter. The 1-body and 2-body terms, which implementation-wise are the only modifications due to the dressed auxiliary construction (v(r) → v (r, q), w(r, r ) → w (rq, r q )) are the same for HF and RDMFT and thus also for dHF and dRDMFT. This means that the validation of dHF that we presented in the S14 previous section at the same time largely validates the implementation of dRDMFT.
Still, we need to analyze and understand the convergence with respect to the basis set in dRDMFT and check for the consistency between the results of RDMFT and dRDMFT in the no-coupling limit.

Basis-set convergence of RDMFT
In RDMFT, we have to perform two minimizations that are interdependent: for the natural orbitals φ i and the natural occupation numbers n i (where always i = 1, . . . , M .) This is done by fixing alternately φ i or n i , while optimizing the other until overall convergence is achieved.
We have the possibility to define different convergence criteria for each minimization, E (which is connected to Λ , see Sec. 2.2) and µ . The latter tests the convergence of the Lagrange multiplier µ that appears in the RDMFT functional to fix the total number of particles in the system. S2 The occupation number optimization routine at iteration step m sets µ = µ m , minimizes the total energy with respect to the n i = n m i and calculates the particle number N m = M i=1 n i . Based on the deviation to the correct system's particle number, µ m+1 is increased or decreased. The routine exits if |µ m − µ m−1 | < µ . For the remainder of this chapter, we set E = µ = 10 −8 .
However, the ES-convergence of RDMFT is not as straightforward as in the HF basis -case.
Contrary to the clear monotonic dependence between ES and the energy that we always found for HF basis , the current RDMFT implementation in Octopus is not strictly variational with respect to the number of basis states. For all tested systems, the energy went down with increasing ES until a certain value ES min and then up again. Therefore, it seems that the interplay of the two minimization processes and the relatively soft types of convergence criteria introduces for high ES such large errors that they exceed the gain of accuracy due to more variational freedom. In Tab. S2, we show the variation of the ground-state energy of He (outside the cavity) for a series ES = {10, 20, . . . , 80}. The energy first decreases until S15 its lowest value for ES min ≈ 40 and then increases again. Table S2: Ground state energies of the He atom (outside the cavity) calculated by RDMFT with parameters mentioned in the main text. The energy goes first down with increasing number of ES until it reaches its lowest value at ES min ≈ 40 and then up again.

ES
Energy ∆E ES = E ES − E ES−10 10 -2.2421837 -20 -2.2426908 −5.1 · 10 −4 30 -2.2427080 −1.7 · 10 −5 40 -2.2427085 −4.9 · 10 −7 50 -2.2427049 +3.6 · 10 −6 60 -2.2427035 +1.3 · 10 −6 70 -2.2426928 +1.1 · 10 −5 80 -2.2426937 −9.9 · 10 −8 This non-strictly variational behavior makes a clear definition of the convergence difficult. However, we find that for every considered system there is an optimal region ES opt = We start with the investigation of the energies and find that ∆E ES,ES < 5 · 10 −6 for all combinations 30 ≤ ES, ES ≤ 60, but ∆E ES,ES > 10 −5 when we choose 30 ≤ ES ≤ 60 and 30 ≤ ES ≤ 60 or 10 ≤ ES ≤ 20. We conclude that the first condition for ES opt is met by the interval 30 ≤ ES ≤ 60 with threshold Eopt = 5 · 10 −6 . For the investigation of the second condition, we depict in Fig. S1 the density deviations ρ ES,ES (x) with ES = 60, the upper boundary of the just found interval and ES = 80, that corresponds to the largest basis set of this example. For a better visibility, the curve for ES = 10 is not shown, but it deviates stronger than all the other curves from both ES . We conclude that ∆ρ ES,ES goes S16 down until ES = 20, independently of the reference. However, for ES ≥ 30 this is not the case any more. We find ∆ρ ES,80 ≈ 10 −4 but ∆ρ ES,60 < 5 · 10 −5 . Additionally, we observe two different types of deviations: The curves ρ ES,80 (x) have a similar form for all 30 ≤ ES ≤ 60, but when we change the reference point to ES = 60, we cannot find pronounced similarities which is what we would expect from fluctuations. When we test also the other possible values for ES , we find ∆ρ ES,ES < 5 · 10 −5 for all 30 ≤ ES, ES ≤ 60. Thus, the second condition for the optimal region is met by the same interval like the first condition with threshold ρopt = 5 · 10 −5 . Therefore, we have ES opt = {ES | 30 ≤ ES ≤ 60} and conclude that the maximum possible accuracy of RDMFT calculations is already reached for ES = 30

{ES | ES
and it is lower than for HF basis . The strong decrease in accuracy of RDMFT in comparison to HF basis suggests that the occupation number minimization adds a significant error to the entire calculation. The cur-S17 rent implementation uses a Broyden-Fletcher-Goldfarb-Shanno algorithm 9 that describes the region around the minimum by a second order Taylor expansion and approximates the corresponding coefficient matrix (Hessian) at every iteration step by a combination of the gradients of the current and the previous iterations, which makes the routine numerically very efficient. It is very difficult to estimate the exact error that is introduced by the method, because of the strong non-linear character of the minimization (especially the interdependence between the φ i and n i optimizations.) For a thorough understanding of this issue, one needs to implement and test different numerical solvers. However, for the purposes of this paper, we consider the current accuracy as sufficient.

Basis-set convergence of dRDMFT and comparison to RDMFT
For the ES-series of dRDMFT, we need to deal with the combination of the inaccuracies introduced by the n i -minimization, explained in the last section and the additional errors due to extra large basis sets that contain many redundant degrees of freedom, that we found for dHF before (Sec. 2.) We know already that a large photon frequency is advantageous in terms of the latter. So let us choose ω = 5.0 to make this error as small as possible.
As expected, it is much harder to define a convergence region in the way we did before for electronic RDMFT. We find the lowest energy at ES = 50, which deviates from ES = 40 and ES = 60 by |E ES=40/60 − E ES=50 | > 10 −5 . So we cannot find a region, where the energy is converged better than 10 −5 . Nevertheless, if we accept an accuracy of ∆E ES < 5 · 10 −5 as sufficient, we find a region as large as 20 ≤ ES ≤ 100 that satisfies the criterion. A look in the density deviations reveals that we can slightly tighten this region and exclude ES = 20, such that we have deviations ∆ρ ES,50 ≈ 10 −4 for all 30 ≤ ES ≤ 100. We found similar values also for other examples and we can conclude that we lose about one order of magnitude of accuracy for dRDMFT in comparison to RDMFT. Careful fine-tuning of the numerical parameters could improve this, but the current accuracy suffices for our purposes.
Finally, we present the consistency check between RDMFT and dRDMFT in the λ = 0 limit. We find |E RDM F T − E dRDM F T | ≈ 10 −5 and max x |ρ RDM F T (x) − ρ dRDM F T (x)| ≈ 10 −4 , which means that the deviations between the levels of theory are of the same order as the maximal accuracy that dRDMFT provides. We conclude that both theories are consistent.
4 Protocol for the convergence of a dHF/dRDMFT calculation We conclude this supporting information with a step-by-step guide for the proper usage of the dressed orbital implementation in Octopus. All the calculations, presented in the main part of this paper were performed according to this protocol 1. Box length L x and spacing dx convergence for the purely electronic part of the system on the level of IP and electronic HF and for the uncoupled photonic system on the level of IP.
• Test the deviations in energy ∆E Lx < E Lx and density ∆ρ Lx < ρ Lx , as explained in Sec. 1. We chose E Lx = 10 −8 and ρ Lx = 10 −5 to exclude any numerical artefacts. However, as the dHF and dRDMFT calculations do typically not reach such precisions, one can relax these criteria in general.
• For the dx-series, test only the deviations in energy ∆E dx < E dx due to the larger density errors. We chose E dx = 10 −8 , but like for the box length, this criterion can be relaxed.
2. Basis size convergence for the HF basis routine with the purely electronic part of the system.
• Compare the converged HF basis and the electronic HF results in energy and density and make sure that both are consistent on their level of accuracy.
3. Basis size convergence of the dressed theory that is wanted (dHF or dRDMFT) in the no-coupling (λ = 0) limit • Perform an ES-series like for HF basis . Note that one needs to expect considerably larger basis sets for the same level of convergence (see Sec. 2.3.2 for details.) • Check consistency of the electronic sub-part of the system with HF basis as mentioned in Sec. 2.3.2. If this check fails drastically, this is very probably due to the violation of the extra exchange symmetry in the photonic coordinates (see Sec. The "fermionization" of matter-photon systems in the main text.) 4. The convergence study is finished with another basis-set convergence for λ > 0. Typically, we also performed another small box length series with the converged basis set to make sure that the coupling does not increase the size of the system crucially such that boundary effects could influence the results.

The bosonic symmetry of the photon wave function
In this appendix we go into a little more detail and show how the mode-representation, which makes the bosonic symmetry explicit, arises and introduce in this setting the usual bosonic density matrices. S4 Instead of starting with the displacement representation we start with the definition of the single-particle Hilbert space and its Hamiltonian. We choose the singleparticle Hilbert space H 1 to consist of M orthogonal states |α . These states are defined by the eigenstates of the Laplacian with fixed boundary conditions and geometry and correspond to the Fourier modes of the electromagentic field. S5,S6 This real-space perspective is a natural choice if one either wants to connect to quantum mechanics and deduce the Maxwell field from gauge independence of the electronic wave function, S5 or when deducing the theory in analogy to the Dirac equation. S7 It is this analogy of Maxwell's equations as a singlephoton wave function with spin 1 that makes the appearance of a bosonic symmetry most explicit when quantizing the theory. S8 We note, however, that in general the concept of a photon wave function can become highly non-trivial. S9 Since we work directly in the dipole approximation we do not go through all the steps of the usual quantization procedure of QED but from the start assume that we have chosen a few of these modes |α (with a certain frequency and polarization) in Coulomb gauge. S6 The single-particle Hamiltonian in this representation is then given byĥ ph = M α=1 ω α |α α| .
Since a total shift of energy does not change the physics and for later reference, we can equivalently useĥ ph = M α=1 ω α + 1 2 |α α|. Therefore, the energy of a single-photon wave function |φ = M α=1 φ(α) |α (corresponding to the classical Maxwell field in Coulomb gauge S8 ) is given by Here we have introduced the single-particle photonic 1RDM γ b (α, β) = φ * (β)φ(α). We can then extend the single-particle space and introduce photonic many-body spaces H N b which are the span of all symmetric tensor products of single-particle states of the form S4,S10 |α 1 , ..., where ℘ goes over all permutations of α 1 , ..., α N b . This construction is completely analogous to the typical construction of the fermionic many-body space with the only difference having minus sings in front of odd permutations. Such a many-body basis is not normalized for bosons, as states can be occupied with more than one particle. Thus, the normalization factor occurs in the corresponding resolution of identity, i.e. 1 = 1 N b ! M α 1 ,...,α N b =1 |α 1 , ..., α N b α 1 , ..., α N b |. This approach is explained in great detail in Ref. S11. An N b -particle Hamiltonian is then given by a sum of individual single-particle Hamiltonians (interactions among the photons will only come about due to the coupling with the electrons.) Introducing for a general N b photon state |φ = 1 √ N b ! M α 1 ,...,α N b =1φ (α 1 , ..., α N b ) |α 1 , ..., α N b withφ(α 1 , ..., α N b ) = 1 √ N b ! α 1 , ..., α N b |φ the corresponding 1RDM according to Eq. (6) as γ b (α, β) = N b α 2 ,...,α N b φ * (β, α 2 , ..., α N b )φ(α, α 2 , ..., α N b ), the energy of that state is given by Such a state can be constructed, for instance, as a permanent of N b single-photon states φ(α). Note further that the 1RDM of an N b photon state obeys N b = α γ b (α, α).
Finally, since we want to have a simplified form of a field theory without fixed number of photons, we make a last step and represent the problem on a Hilbert space with indetermined number of particles, i.e., a Fock space. By defining the vacuum state |0 , which spans the one-dimensional zero-photon space, the Fock space is defined by a direct sum of N b -photon Introducing the ladder operators between the different photonnumber sectors of F by S4 a + α |α 1 , ..., α N b = |α 1 , ..., α N b , α â α |α 1 , ..., α N b = N b k=1 δ α k ,α |α 1 , ..., α k−1 , α k+1 , ..., α N b with the usual commutation relations, we can lift the single-particle Hamiltonian to the full Fock space and arrive at Eq. (3). The Fock space 1RDM for a general Fock space wave