CERN-LHC. Search for the pair production of photon-jets---collimated groupings of photons---with the ATLAS detector. Highly collimated photon-jets can arise from the decay of new, highly boosted particles that can decay to multiple photons collimated enought to be identified in the electromagnetic calorimeter as a single, photon-like energy cluster. Data from proton–proton collisions at a centre-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 36.7 fb$^{-1}$, were collected in 2015 and 2016. Upper limits are placed on $\sigma\times \mathcal{B}(X\rightarrow aa)\times \mathcal{B}(a\rightarrow\gamma\gamma)^2$ and $\sigma\times \mathcal{B}(X\rightarrow aa)\times \mathcal{B}(a\rightarrow 3\pi^0)^2$ for 200 GeV < $m_X$ < 2TeV and $m_a$ < 10 GeV. Tables 8 to 35 are provided to allow the recasting of the cross-section upper limits to different signal models predicting final states with photon-jets. These tables present the selection efficiency (before categorisation) $\varepsilon_{\gamma_R}(E_\mathrm{T},\eta)$ for reconstructed photons originating from a photon-jet, and the fraction $f_{\gamma_R}(E_\mathrm{T},\eta)$ of reconstructed photons with a value of the shower shape variable $\Delta E$ lower than the threshold. The fiducial region is defined as: - $E_\mathrm{T,1}>0.4\times m_X$ - $E_\mathrm{T,2}>0.3\times m_X$ - $|\eta_i| < 2.37 (i=1,2)$ (excluding $1.37 < |\eta_i| <1.52$) where $E_\mathrm{T,1}, \eta_1$ ($E_\mathrm{T,2}, \eta_2$) are the transverse energy and the pseudorapidity of the $a$ particle with the higher (the lower) transverse energy, respectively. For a resonance particle $X$ decaying into a pair of photon-jets via $X\rightarrow aa$, the total selection efficiency, $\varepsilon$, and the fraction of events in the low-$\Delta E$ category, $f$, can be computed by integrating over the p.d.f. of $(E_\mathrm{T,1},\eta_1,E_\mathrm{T,2},\eta_2)$ with the following procedure: - apply the fiducial cuts to the two $a$ particles - compute $\varepsilon$ from the integration of $\varepsilon_{\gamma_R}(E_\mathrm{T,1},\eta_1) \cdot \varepsilon_{\gamma_R}(E_\mathrm{T,2},\eta_2)$ - compute $f$ from the integration of $\varepsilon_{\gamma_R}(E_\mathrm{T,1},\eta_1) \cdot \varepsilon_{\gamma_R}(E_\mathrm{T,2},\eta_2) \cdot f_{\gamma_R}(E_\mathrm{T,1},\eta_1) \cdot f_{\gamma_R}(E_\mathrm{T,2},\eta_2)$ divided by $\varepsilon$ With the resulting value of $f$ for a given value of $m_X$, the 95% CL observed upper limit on the visible cross-section (i.e. $\sigma\times \mathcal{B}\times\varepsilon$) can be taken from Table 7, which is considered to be model-independent. The corresponding upper limit on the cross-section times branching ratios, $\sigma \times \mathcal{B}$, can be computed by dividing the obtained visible cross-section by $\varepsilon$. The estimation procedure described above is validated by comparing the results for the benchmark signal scenario decaying via $X\rightarrow aa\rightarrow 4\gamma$ with the results presented in the paper (i.e. Table 3). It is found that the two results agree within 20%, and the result with the estimation procedure described above gives lower values. The main difference is found for large values of the mass ratio, $0.005<m_a/m_X$. This is because, for larger values of $m_a/m_X$, the width of the distribution of the reconstructed diphoton mass $m_{\gamma_{R}\gamma_{R}}$ increases. This is caused by the wider angular separation between the photons inside a photon-jet for larger $m_a/m_X$, leading to a greater part of the energy of the shower leaking out of the window defined in the cells of the electromagnetic calorimeter to collect energy for the photon reconstruction.