# 3D Fermi surfaces from charge order in layered CsV3Sb5

Xiangwei Huang  Laboratory of Quantum Materials (QMAT), Institute of Materials
(IMX),École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland  Chunyu Guo∗ Laboratory of Quantum Materials (QMAT), Institute of
Materials (IMX),École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Max Planck Institute for the Structure and Dynamics of
Matter, 22761 Hamburg, Germany  Carsten Putzke  Laboratory of Quantum
Materials (QMAT), Institute of Materials (IMX),École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Max Planck Institute for the
Structure and Dynamics of Matter, 22761 Hamburg, Germany  Martin Gutierrez-
Amigo  Department of Physics, University of the Basque Country (UPV/EHU),
48080 Bilbao, Spain  Centro de Física de Materiales (CSIC-UPV/EHU), 20018
Donostia-San Sebastian, Spain  Yan Sun  Max Planck Institute for Chemical
Physics of Solids, 01187 Dresden, Germany  Maia G. Vergniory  Max Planck
Institute for Chemical Physics of Solids, 01187 Dresden, Germany  Donostia
International Physics Center, 20018 Donostia-San Sebastian, Spain  Ion Errea
Donostia International Physics Center, 20018 Donostia-San Sebastian, Spain
Centro de Física de Materiales (CSIC-UPV/EHU), 20018 Donostia-San Sebastian,
Spain  Fisika Aplikatua Saila, Gipuzkoako Ingeniaritza Eskola, University of
the Basque Country (UPV/EHU), 20018 Donostia-San Sebastian, Spain  Dong Chen
Max Planck Institute for Chemical Physics of Solids, 01187 Dresden, Germany
Claudia Felser  Max Planck Institute for Chemical Physics of Solids, 01187
Dresden, Germany  Philip J. W. Moll† Laboratory of Quantum Materials (QMAT),
Institute of Materials (IMX),École Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland  Max Planck Institute for the Structure and
Dynamics of Matter, 22761 Hamburg, Germany

## Abstract

The cascade of electronic phases in CsV3Sb5 raises the prospect to disentangle
their mutual interactions in a clean, strongly interacting Kagome lattice.
When the Kagome planes are stacked into a crystal, its electronic
dimensionality encodes how much of the Kagome physics and its topological
aspects survive. The layered structure of CsV3Sb5 reflects in Brillouin-zone-
sized quasi-2D Fermi surfaces and a significant transport anisotropy. Yet here
we demonstrate that CsV3Sb5 is a three-dimensional metal within the charge-
density-wave (CDW) state. Small 3D pockets play a crucial role in its low-
temperature magneto- and quantum transport. Their emergence at $T_{CDW}\sim
93K$ results in an anomalous sudden increase of the in-plane magnetoresistance
by 4 orders of magnitude. The presence of these 3D pockets is further
confirmed by quantum oscillations under in-plane magnetic fields -
demonstrating their closed nature. These results emphasize the impact of
interlayer coupling on the Kagome physics in 3D materials.

## introduction

When structurally layered materials host strong electronic correlations, their
effective electronic dimensionality is key to understanding their microscopic
physics. The anisotropic Fermi surface then sets the canvas on which
correlated ground states unfold, such as magnetism, charge-density-waves (CDW)
or superconductivity. Cuprate and pnictide high-temperature superconductors
are prime examples of layered materials in which reduced dimensionality
defines the superconducting stateBednorz and Müller (1986); Kamihara  _et al._
(2006). The anisotropy is further a critical parameter to describe the vortex
formation and the orbital limit of the upper critical field. Recently, metals
hosting planes with Kagome nets have attracted significant attention due to
their topological phases as well as potential correlation effectsKida  _et
al._ (2011); Kang  _et al._ (2020); Ortiz  _et al._ (2019). The latter have
been argued to arise from small-bandwidth bands (flat bands) that are
predicted in phenotypical two-dimensional Kagome modelsO’Brien  _et al._
(2010). Flat bands naturally enhance correlation effects when the Coulomb
energy exceeds the kinetic energy and they lead to a divergent density of
states, an entropic catastrophe which is commonly avoided by the formation of
electronically ordered states at low temperatures. However when the Kagome
nets are stacked into a three-dimensional crystal, the existence of flat bands
depends on the strength of the interlayer hybridization and its impact on the
band dispersion. Therefore in order to explore how much Kagome physics
survives in a 3D metal, it is necessary to investigate the electronic
structure and especially its effective dimensionality.

Here we investigate the electronic dimensionality of CsV3Sb5 [Fig. 1(a) and
(b)], a stacked Kagome crystal, and its influence on the rich physics it hosts
such as the non-trivial band topology, superconductivity and CDW order. The
non-magnetic Kagome nets formed by Vanadium atoms result in a symmetry-
enforced pair of $\mathbb{Z}_{2}$ topological bands with opposite Chern
numbers Ortiz  _et al._ (2019, 2020); Fu  _et al._ (2021); Chen  _et al._
(2022). The hybridization with Sb atoms effectively expands the band width and
therefore the flat band gains sizable dispersion yet it leaves multiple bulk
Dirac points intact, as observed by angle-resolved photoemission spectroscopy
(ARPES)Ortiz  _et al._ (2020); Luo  _et al._ (2021). Upon cooling, the Kagome
lattice undergoes a CDW transition at $T_{CDW}$ = 93 K Yu  _et al._ (2021a);
Mu  _et al._ (2021); Liang  _et al._ (2021). The 2$\times$2 reconstruction
within the Kagome plane is accompanied by a $\pi$ phase shift of the CDW
across an atomic step edge as observed by STM experiments Zhao  _et al._
(2021); Liang  _et al._ (2021), suggesting a band reconstruction along the
$c$-direction. X-ray studies further suggest a 2$\times$2$\times$4
superlattice formation Wu  _et al._ (2022); Stahl  _et al._ (2022) with long
range out-of-plane coherence, again highlighting the importance of out-of-
plane coupling in this system. At even lower temperatures of $T_{c}\sim 2.8$
K, superconductivity appearsOrtiz  _et al._ (2020). The large upper critical
field anisotropy ($\sim$ 9) is a result of its anisotropic superconducting
propertiesNi  _et al._ (2021). Meanwhile the non-trivial band topology may
give rise to topological superconductivity. Indeed, STM has uncovered a zero-
bias peak in the vortex cores which, while not conclusive, is compatible with
a scenario of bound Majorana statesLiang  _et al._ (2021).

## Results

The single-particle band structure from ab-initio calculations serves as a
starting point to investigate the dimensionality of this strongly correlated
compound. The electronic structure of CsV3Sb5 in the high-symmetry state at
room temperature is calculated using density functional theory (DFT) methods,
the computational details can be found in the Supplementary materialsGiannozzi
_et al._ (2017); Perdew  _et al._ (1996); Dal Corso (2014); Methfessel and
Paxton (1989); Mostofi  _et al._ . Multiple bands are found at the Fermi level
[Fig. 1(c)], consistent with the previous reportsOrtiz  _et al._ (2019, 2020);
Fu  _et al._ (2021). Accordingly, multiple Fermi surfaces of cylindrical shape
are expected in the Brillouin zone [Fig. 1(d)]. While the layered structure is
reflected in the electron-like cylinders centered at $\Gamma$ and K, their
sizable dispersion signals pronounced interlayer coupling. In absence of spin-
orbit coupling (SOC), the cylinder at the Brillouin zone boundary (H to K)
would surround a nodal line, yet SOC breaks this degeneracy gapping the nodal
line.

The band structure of CsV3Sb5 indicates a clear out-of-plane dispersion,
giving rise to warping of the cylindrical Fermi surfaces. Experimentally the
electronic dimensionality can be explored via the resistivity anisotropy. As
common in structurally layered materials, CsV3Sb5 crystals grow as thin
platelets along the Kagome plane Ortiz  _et al._ (2019). While this crystal
morphology lends itself to in-plane resistivity measurements, quantitative
out-of-plane transport poses a well-known challenge. Focused Ion Beam (FIB)
milling can be used to prepare micron-sized $c$-direction bars with well-
defined geometries to quantify the resistivityMoll (2018). CsV3Sb5 like other
soft CDW-compounds is susceptible to mechanical strainSong  _et al._ (2021);
Du  _et al._ (2021), we suspend the microstructured sample (purple) in free
space. It is mechanically and electrically connected to a supportive silicon
frame only via thin, gold-coated SiNx membrane microsprings (gold) [Fig. 1(e)
and (f)]. These structures exert only minimal residual pressure on the sample
($\approx$ 9.8 bar, see supplement), thus allowing an intrinsic evolution of
the ordered states. Moreover, a low-voltage polish step at 5 kV reduces the
amorphization layer into the nm-rangeKelley  _et al._ (2013) which minimizes
its possible influence.

In these suspended structures, the in-plane resistivity decreases from 80
$\mu\Omega$cm at 300 K to 0.3 $\mu\Omega$cm at $T_{c}=2.8$ K. The high
residual resistivity ratio (RRR) of 250 is comparable to bulk samplesChen  _et
al._ (2022), evidencing the high quality of the microstructure fabrication
(Fig. 2). The in-plane resistivity displays a small yet sharply defined
discontinuity at the CDW transition at $T_{CDW}$ which is consistent with
previous reportsChen  _et al._ (2022); Ortiz  _et al._ (2019); Xiang  _et al._
(2021). The out-of-plane resistivity, $\rho_{c}$, is comparably larger yet
shows an overall similar metallic temperature dependence. It decreases from
$\rho_{c}$(300K) $\sim$ 1.9 m$\Omega$cm to $\rho_{c}$($T_{c}$) $\sim$ 33
$\mu\Omega$cm, corresponding to a lower RRR of 60. These results are
qualitatively reflected in the challenging measurements of $\rho_{c}$ on bulk
crystals Ortiz  _et al._ (2020); Xiang  _et al._ (2021), yet deviate
quantitatively. In comparison to the in-plane resistivity, $\rho_{c}$ shows an
upwards jump at $T_{CDW}$ [Fig. 2(c)]. This jump
$(\rho_{c}(T_{CDW}^{+})-\rho_{c}(T_{CDW}^{-}))/\rho_{c}(T_{CDW}^{+})\approx
12\%$ is much more pronounced compared to the change of $\approx-3\%$ for in-
plane resistivity. This relative difference again emphasizes the importance of
the CDW in the out-of-plane direction.

Further information about the electronic dimensionality is contained in the
temperature dependence of the anisotropy. At high temperatures above
$T_{CDW}$, the picture of an anisotropic metal emerges
[$\rho_{c}/\rho_{ab}(T=300~{}K)\approx 24$], which falls close to the Fermi
velocity anisotropy when averaged over the entire Fermi surface
($|v_{F,\parallel}|/|v_{F,\perp}|\sim 20$, [Fig. 2(d)]). The dominant carrier
density originates from the Brillouin-zone-sized hexagonal Fermi surfaces with
weak $k_{z}$ dispersion. These feature six almost flat surfaces parallel to
the crystalline $a$-directions [Fig. 1(b)], which play an important role in
the magneto-transport as discussed later.

Further lowering the temperature beyond the jump at $T_{CDW}$ increases the
anisotropy which eventually saturates around 100 at $T_{c}$, suggesting
dominant in-plane electronic transport at low temperatures. Yet $\rho_{c}$
reaches a low value of 33 $\mu\Omega cm$ at $T_{c}$, which signals metallic
three-dimensional transport. The first-order derivative of the resistivity
anisotropy reveals further, more subtle changes in the material [Fig. 2(e)].
Besides the clear spike at $T_{CDW}$, a local minimum occurs at around 70 K.
This temperature coincides with the onset of anisotropy in the muon spin
depolarization rate that has been associated with time-reversal-symmetry
breaking (TRSB)Yu  _et al._ (2021b). The TRSB state has been theoretically
proposed to arise from an effective orbital current loop fluxDenner  _et al._
(2021); Feng  _et al._ (2021), and a modification of a magnetic scattering
channel would be a natural connection between this experiment and ours.

This further emphasizes the emerging question about its effective electronic
dimensionality. The large resistivity anisotropy is compatible with an
effective 2D description, which may be captured by a simplified model based on
a 2D Kagome lattice. Yet to discuss the effective dimensionality, one has to
carefully distinguish between the transport dimensionality and the
dimensionality of the Fermi surface. In general, the transport anisotropy is
determined by the Fermi velocity $v_{F}(k)$ and the scattering time $\tau(k)$
distribution while the Fermi surface topology encodes another aspect of
electronic dimensionality. For example, a cylindrical Fermi surface with
strong warping naturally features low transport anisotropy. Here we
demonstrate that CsV3Sb5 is a 3D Kagome metal as the emergent small, closed
Fermi pockets play a crucial role in its magneto- and quantum transport
properties.

The magnetoresistance provides further information about the electronic
dimensionality (Fig. 3). Most intriguing is the complete absence of transverse
magnetoresistance for orthogonal in-plane fields and currents at any
temperature above $T_{CDW}$. The field-independent noise floor provides an
upper bound for the magnetoresistance, at ($\rho_{a}$(18 T)-$\rho_{a}$(0
T))/$\rho_{a}$(0 T) $<$ 10-5 at 120 K. Even rather unsuspicious metals such as
Cu commonly show transverse magnetoresistance around $1\%$ (20T) at room
temperature and $100\%$ (20T) at 100KDe Launay  _et al._ (1959). In contrast,
for all other orthogonal current and field orientations in CsV3Sb5, the
magnetoresistance increases with decreasing temperature following a
conventional semiclassical scaling of the magneto-transport by
$\omega_{c}\tau$. Given the complex Fermi surface and putative incipient TRSB
order, such absent magnetoresistance for one particular direction appears
quite exotic.

The reconstruction of the main cylindrical Fermi surfaces at $T_{CDW}$ offers
a natural explanation for such behavior. When a current is applied in the
plane, the conductivity is dominated by the weakly warped main Fermi surfaces
[Fig. 3(d)]. Such flatness ensures that the quasiparticle velocity remains
approximately unchanged by the Lorentz force. This behavior is commonly
observed in quasi-2D materials when only cylindrical Fermi surfaces are
presentKorin-Hamzić  _et al._ (2002); Wosnitza (2006). For magnetic field or
current applied along any other direction, magnetoresistance is bound to
appear as the Lorentz force significantly impacts the carrier trajectories.
These results consistently demonstrate that CsV3Sb5 is electronically quasi-2D
above $T_{CDW}$.

However, this conclusion no longer holds in the charge ordered phase. Though
the exact type of CDW reconstruction is still under debate, all different
scenarios predict the emergence of small closed 3D pockets as a direct
consequence of zone foldingOrtiz  _et al._ (2021). Compared to the quasi-2D
Fermi surfaces, these small pockets support closed cyclotron orbits which
leads to the sudden recovery of magnetoresistance [Fig. 3(d)]. The direct
connection between the abrupt change of magnetoresistance and Fermi surface
reconstruction further demonstrates the importance of the small, closed
pockets to the 3D electronic transport properties of CsV3Sb5, especially at
low temperatures.

Such small, 3D-like pockets at low temperature should be visible in quantum
oscillations. Four low frequency quantum oscillations have been reported for
out-of-plane fields and associated with these small 3D pocketsOrtiz  _et al._
(2021); Yu  _et al._ (2021a). However, these quantum oscillations were
observed to quickly disappear as the field is rotated towards the plane, which
precluded the distinction between a closed, 3D pocket and an open, 2D
cylinder. The weakly perturbed, suspended microbars now allow to track this
frequency over the entire angle range, firmly establishing them as closed 3D
pockets within the density wave phase (Fig. 4). At low angles, the oscillation
frequencies observed in our structures are consistent with these previous
reports. A likely scenario for the loss of amplitude in conventional
experiments involves stacking disorder along $c$-direction due to residual
strain, which reduces the quantum-coherent transport along the out-of-plane
direction. This is consistent with the larger RRR and lower resistivity
observed in the microstructures compared to macroscopic crystalsOrtiz  _et
al._ (2020); Xiang  _et al._ (2021). The angular dependence of the oscillation
frequencies is well described by a simple ellipsoidal model $F_{i}\propto
1/\sqrt{sin^{2}(\theta)+(1+A_{i}^{2})cos^{2}(\theta)}$. The factor $A_{i}$
characterizes the anisotropy of the i-th Fermi surface [Fig. 4(d)]. These
results further support the presence of 3D pockets as suggested by the unusual
magneto-transport. The emergence of 3D Fermi pockets due to CDW reconstruction
has already been observed in other CDW materials Knowles  _et al._ (2020); Tam
_et al._ (2022) as a 3D zone-folding of Brillouin-zone-sized Fermi surface
naturally results in small Fermi surface pockets. These pockets are connected
by the 3/4 Bragg wavevector in the 3Q directions which results in an exotic
pairing density wave phase at low temperatureZhou and Wang (2021), as observed
in STM experiments Chen  _et al._ (2021).

In conclusion, the electronic anisotropy and Fermiology of charge-ordered
CsV3Sb5 are characterized by the coexistence of a high-mobility in-plane
electronic system with small 3D pockets. They emerge when the charge order is
established, demonstrating an electronic dimensionality crossover at
$T_{CDW}$. Therefore the interlayer coupling must be taken into account for
theoretical modelling as considerations based on 2D Kagome lattices may miss
important aspects of this material. It will be highly interesting to explore
how the cascades of correlated electronic states are influenced by the
coupling between Kagome nets.

## References

## References

  * Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Z. Phys. B Condens. 64, 189 (1986).
  * Kamihara  _et al._ (2006) Y. Kamihara, H. Hiramatsu, M. Hirano, R. Kawamura, H. Yanagi, T. Kamiya, and H. Hosono, J. Am. Chem. Soc. 128, 10012 (2006).
  * Kida  _et al._ (2011) T. Kida, L. Fenner, A. Dee, I. Terasaki, M. Hagiwara, and A. Wills, J. Condens. Mat. Phys. 23, 112205 (2011).
  * Kang  _et al._ (2020) M. Kang, S. Fang, L. Ye, H. C. Po, J. Denlinger, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, J. G. Checkelsky, _et al._ , Nat. Commun. 11, 4004 (2020).
  * Ortiz  _et al._ (2019) B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, _et al._ , Phys. Rev. Mater. 3, 094407 (2019).
  * O’Brien  _et al._ (2010) A. O’Brien, F. Pollmann, and P. Fulde, Phys. Rev. B 81, 235115 (2010).
  * Ortiz  _et al._ (2020) B. R. Ortiz, S. M. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, _et al._ , Phys. Rev. Lett. 125, 247002 (2020).
  * Fu  _et al._ (2021) Y. Fu, N. Zhao, Z. Chen, Q. Yin, Z. Tu, C. Gong, C. Xi, X. Zhu, Y. Sun, K. Liu, and H. Lei, Phys. Rev. Lett. 127, 207002 (2021).
  * Chen  _et al._ (2022) D. Chen, B. He, M. Yao, Y. Pan, H. Lin, W. Schnelle, Y. Sun, J. Gooth, L. Taillefer, and C. Felser, Phys. Rev. B 105, L201109 (2022).
  * Luo  _et al._ (2021) Y. Luo, S. Peng, S. M. Teicher, L. Huai, Y. Hu, B. R. Ortiz, Z. Wei, J. Shen, Z. Ou, B. Wang, _et al._ , arXiv:2106.01248 (2021).
  * Yu  _et al._ (2021a) F. Yu, T. Wu, Z. Wang, B. Lei, W. Zhuo, J. Ying, and X. Chen, Phys. Rev. B 104, L041103 (2021a).
  * Mu  _et al._ (2021) C. Mu, Q. Yin, Z. Tu, C. Gong, H. Lei, Z. Li, and J. Luo, Chin. Phys. Lett. 38, 077402 (2021).
  * Liang  _et al._ (2021) Z. Liang, X. Hou, F. Zhang, W. Ma, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, _et al._ , Phys. Rev. X 11, 031026 (2021).
  * Zhao  _et al._ (2021) H. Zhao, H. Li, B. R. Ortiz, S. M. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Nature 599, 216 (2021).
  * Wu  _et al._ (2022) S. Wu, B. R. Ortiz, H. Tan, S. D. Wilson, B. Yan, T. Birol, and G. Blumberg, Phys. Rev. B 105, 155106 (2022).
  * Stahl  _et al._ (2022) Q. Stahl, D. Chen, T. Ritschel, C. Shekhar, E. Sadrollahi, M. C. Rahn, O. Ivashko, M. v. Zimmermann, C. Felser, and J. Geck, Phys. Rev. B 105, 195136 (2022).
  * Ni  _et al._ (2021) S. Ni, S. Ma, Y. Zhang, J. Yuan, H. Yang, Z. Lu, N. Wang, J. Sun, Z. Zhao, D. Li, _et al._ , Chin. Phys. Lett. 38, 057403 (2021).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys. Condens. Mat. 29, 465901 (2017).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Dal Corso (2014) A. Dal Corso, Comp. Mat. Sci. 95, 337 (2014).
  * Methfessel and Paxton (1989) M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989).
  * (22) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Moll (2018) P. J. W. Moll, Annu. Rev. Condens. Mat. Phys. 9, 147 (2018).
  * Song  _et al._ (2021) B. Song, X. Kong, W. Xia, Q. Yin, C. Tu, C. Zhao, D. Dai, K. Meng, Z. Tao, Z. Tu, _et al._ , arXiv:2105.09248 (2021).
  * Du  _et al._ (2021) F. Du, S. Luo, B. R. Ortiz, Y. Chen, W. Duan, D. Zhang, X. Lu, S. D. Wilson, Y. Song, and H. Yuan, Phys. Rev. B 103, L220504 (2021).
  * Kelley  _et al._ (2013) R. Kelley, K. Song, B. Van Leer, D. Wall, and L. Kwakman, Microsc. Microanal. 19, 862 (2013).
  * Xiang  _et al._ (2021) Y. Xiang, Q. Li, Y. Li, W. Xie, H. Yang, Z. Wang, Y. Yao, and H.-H. Wen, Nat. Commun. 12, 6727 (2021).
  * Yu  _et al._ (2021b) L. Yu, C. Wang, Y. Zhang, M. Sander, S. Ni, Z. Lu, S. Ma, Z. Wang, Z. Zhao, H. Chen, _et al._ , arXiv:2107.10714 (2021b).
  * Denner  _et al._ (2021) M. M. Denner, R. Thomale, and T. Neupert, Phys. Rev. Lett. 127, 217601 (2021).
  * Feng  _et al._ (2021) X. Feng, K. Jiang, Z. Wang, and J. Hu, Sci. Bull. 66, 1384 (2021).
  * De Launay  _et al._ (1959) J. De Launay, R. Dolecek, and R. Webber, J. Phys. Chem. Solids 11, 37 (1959).
  * Korin-Hamzić  _et al._ (2002) B. Korin-Hamzić, M. Basletić, and K. Maki, Europhysics Letters (EPL) 59, 298 (2002).
  * Wosnitza (2006) J. Wosnitza,  _Fermi surfaces of low-dimensional organic metals and superconductors_ , Vol. 134 (Springer, 2006).
  * Ortiz  _et al._ (2021) B. R. Ortiz, S. M. Teicher, L. Kautzsch, P. M. Sarte, J. P. Ruff, R. Seshadri, and S. D. Wilson, Phys. Rev. X 11, 041030 (2021).
  * Knowles  _et al._ (2020) P. Knowles, B. Yang, T. Muramatsu, O. Moulding, J. Buhot, C. J. Sayers, E. Da Como, and S. Friedemann, Phys. Rev. Lett. 124, 167602 (2020).
  * Tam  _et al._ (2022) C. C. Tam, M. Zhu, J. Ayres, K. Kummer, F. Yakhou-Harris, J. R. Cooper, A. Carrington, and S. M. Hayden, Nature Communications 13, 570 (2022).
  * Zhou and Wang (2021) S. Zhou and Z. Wang, arXiv:2110.06266 (2021).
  * Chen  _et al._ (2021) H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, S. Ma, S. Ni, H. Zhang, Q. Yin, C. Gong, Z. Tu, H. Lei, H. Tan, S. Zhou, C. Shen, X. Dong, B. Yan, Z. Wang, and H.-J. Gao, Nature 599, 222 (2021).

Acknowledgements

Funding:  This work was funded by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (MiTopMat
- grant agreement No. 715730). This project received funding by the Swiss
National Science Foundation (Grants No. PP00P2_176789). M.G.V., I. E. and
M.G.A. acknowledge the Spanish Ministerio de Ciencia e Innovacion (grant
PID2019-109905GB-C21). M.G.V. thanks support to Programa Red Guipuzcoana de
Ciencia Tecnología e Innovación 2021 No. 2021-CIEN-000070-01 Gipuzkoa Next and
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) GA
3314/1-1 – FOR 5249 (QUAST). This work has been supported in part by Basque
Government grant IT979-16. This work was also supported by the European
Research Council Advanced Grant (No. 742068) “TOPMAT”, the Deutsche
Forschungsgemeinschaft (Project-ID No. 258499086) “SFB 1143”, and the DFG
through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology
in Quantum Matter ct.qmat (EXC 2147, Project-ID No. 39085490).

Competing Interests The authors declare that they have no competing financial
interests.

Figure 1: Side(a) and top(b) view of crystal structure of CsV3Sb5. It consists
with layers of Cs, V and Sb atoms which form the Kagome lattice plane. (c)
Electronic structure of CsV3Sb5 calculated by density-functional theory. There
exist multiple electronic bands across the Fermi level. (d) Fermi surfaces in
the Brillouin zone, including several types of Fermi surfaces and their
symmetric copies. (e) Illustration of membrane-mount CsV3Sb5 device. The
sample is suspended with soft SiNx springs, which ensures the minimum strain
effect. (f) Scanning-electron-microscope (SEM) image of the Membrane device.
The 250 by 250 $\mu$m membrane window was patterned into several branches of
springs with FIB. Figure 2: (a) Temperature dependence of the out-of-
plane($\rho_{c}$) and in-plane resistivity($\rho_{ab}$). (b) An abrupt jump
occurs in the temperature dependence of $\rho_{c}$, which corresponds to the
charge-density-wave transition temperature $T_{CDW}$. Consistently $\rho_{ab}$
also displays a weak discontinuity. (c) Enlarged view of (a) within the low
temperature range which demonstrates a clear superconducting transition at
$T_{c}$ = 2.8 K. (d) Temperature dependence of resistivity anisotropy
($\rho_{c}/\rho_{ab}$). The dashed lines stand for the averaged Fermi velocity
anisotropy for all Fermi surfaces. (e) First-order derivative of temperature
dependence of resistivity anisotropy [d($\rho_{c}/\rho_{ab}$)/d$T$]. At the
charge-density-wave transition temperature ($T_{CDW}$) a sharp peak can be
observed, as well as a broad hump at $T$” $\approx$ 14 K.

Figure 3: Field dependence of out-of-plane (a) and in-plane (b)
magnetoresistance ratio [MR = ($\rho$ (18 T)-$\rho$ (0 T))/$\rho$ (0 T)] with
field applied along a and c(ab)-axis. (c) Temperature dependence of
magnetoresistance ratio at $B$ = 18 T with current along ab-axis and field
along a-axis which are defined in the inset. (d) Illustration of Fermi surface
reconstruction due to charge ordering. The abrupt appearance of
magnetoresistance below $T_{CDW}$ is a direct consequence of the 3D pockets
appears only after the CDW reconstruction.

Figure 4: (a) Fast-Fourier-transformation(FFT) spectrum with different field
directions for Shubnikov-de-Hass(SdH) oscillations measured with current
applied along c-axis. Here $\theta$ stands for the angle between a-axis and
field direction. (b) SdH oscillations ($\rho_{osc}$) at $\theta$ = 6, 0 and
-10∘. Here $\rho_{osc}=\Delta\rho/\rho_{BG}$, with $\Delta\rho$ the
oscillatory part of the magneto-resistivity, and $\rho_{BG}$ a polynomial fit
to the magneto-resistivity background. (c) Corresponding FFT spectrum with the
star symbols indicate the identified peaks. (d) Angular dependence of SdH
oscillation frequencies. The error bar is defined as full width at half
maximum of the peaks in the FFT spectrum.


# Light Scattering on Two Parallel Subwavelength Cylinders

Alex Verevkin  Institute of Automation and Electrometry, Siberian Branch,
Russian Academy of Sciences, 1 Koptjug Avenue, Novosibirsk 630090, Russia
Alex Bereza  Institute of Automation and Electrometry, Siberian Branch,
Russian Academy of Sciences, 1 Koptjug Avenue, Novosibirsk 630090, Russia
David Shapiro  Institute of Automation and Electrometry, Siberian Branch,
Russian Academy of Sciences, 1 Koptjug Avenue, Novosibirsk 630090, Russia

###### Abstract

The plasmon resonance has found important application in various systems,
e.g., nanoantennas, solar panels, refractive index sensors. Unfortunately, a
few analytical solutions for such systems are known. The work aims to find a
solution for scattering by a plane electromagnetic wave on two parallel
cylinders. When their diameters and the gap between them are less than the
radiation wavelength, the quasistatic approximation is valid. We build up a
conformal transformation that maps Cartesian into bipolar orthogonal
coordinates and represent the scattered field as a decomposition by
eigenfunctions of the Laplace operator. The near-field intensity distribution
is shown to coincide with numerical calculation performed by COMSOL
Multiphysics. Green’s function obtained allows one to find also higher
decomposition orders. Compared to previous studies, we treat the general case
of cylinders with different diameters and dielectric constants.

##  I Introduction

One of the most promising and exciting effects in modern nano-optics
girard2005 ; shalaev2006  is a plasmon resonance, the excitation of electron
plasma oscillations in a metal under the electromagnetic waves. These
oscillations lead to the local field enhancement when the incident field could
increase by several orders stockman2011 ; novotny2012 . Moreover, the light
localization at the nanoscale can be achieved by an optical resonance in pure
dielectric structures with no metallic elements devilez2015 ; krasnok2018 .
Therefore, it is possible to obtain a considerable intensity in a particular
geometrical structure.

The applications of plasmon resonances are various. For example, plasmonic
structures in the photovoltaic panel can significantly enhance its efficiency
by increasing the energy conversion ratio of the incident wave. Usually, one
increases the intensity by confining the semiconductors’ energy, where plasmon
structures create additional electron-hole pairs enrichi2018 . The consequent
application is the optical sensors of the refraction index based on the
plasmon resonance effect homola2008 ; shalabney2011 ; piliarik2012 . By
analyzing the spectrum of reflection and shift of the resonance peaks in the
angular dependence, one can reach high values of the resolution and
sensitivity xu2019  in conjunction with their comparatively small size and the
possibility of repeated use. For instance, one improves Bragg sensors covering
the fiber surface by nanoparticles bialiayeu2012 ; arasu2016 . However, all
these applications face a significant obstacle. Due to the non-trivial
geometry, it is almost impossible to find analytical solutions for a scattered
amplitude. Current methods for this problem analysis are numerical: Finite-
Difference Time-Domain method (FDTD) yee1966numerical , Discrete Dipole
Approximation (DDA) purcell73 ; Draine08 , Boundary Element Method (BEM)
brebbia2016boundary , Finite Element Method (FEM) zienkiewicz2005finite .
However, each method has limitations and numerical errors, which can not be
quickly evaluated and determined. Therefore, new analytical solutions must
assist in choosing an optimal configuration for the plasmon structure. At
least it is helpful to verify the numerical simulations.

The present paper study the scattering of plane waves by two parallel
cylinders. The problem is two-dimensional, and the aim is to find the field
distribution using a conformal transformation. Such a quasistatic approach is
valid for subwavelength particles VKS7 . A conformal transformation for two
cylinders, the bipolar coordinates, had been found earlier for the symmetric
situation of equal cylinders vorobev2010 . Below, we expand the method to
cylinders with different diameters and dielectric constants. Ref. lei2010  has
considered the system of two wires for a specific case, when their cross-
sections have a mutual point. The authors call it the ”kissing nanowires.”
Energy concentration and confinement occur in the singularity, a point of
contact. Surface plasmon modes propagate toward the singularity, place of the
energy accumulation. There is an alternative way to obtain analytical
approximate solutions, the Born approximation modified for dielectric
particles bereza2017 ; bereza2019 ; PhysRevRes2020 ; PhysRevA.104.063514 ;
Ustimenko2022 . Chains of parallel infinite cylinders also permits approximate
analytical solutions in the form of decomposition over cylindrical waves
belan2015 ; lee2016 .

We introduce the geometry in Sect. II and discuss the conformal transformation
in Sect. III. Next, in Sect. IV, we find the zero-order solution in the form
of decomposition by eigenfunctions. The near-field intensity distribution is
compared with numeical calculation. Sect. V summarizes the conclusions. We
present Green’s function, allowing one to calculate the next orders, in
Appendix.

##  II Laplace equation

Fig. 1 shows the cross-section of considered geometry. Parallel cylinders of
radius $R_{1}$ and $R_{2}$, with dielectric permittivity $\varepsilon_{1}$ and
$\varepsilon_{2}$, respectively, interact with plane electromagnetic wave of
amplitude $\vec{E}_{0}$ and wavevector $\vec{k}$. A gap between two cylinders
has width $\Delta$. Axis $y$ of the system goes through centers of circles,
while $z$ is parallel to the cylinder axes. Validity condition is that the
wavelength of incident wave $\lambda$ should be greater than each
characteristic dimension of our system

$kb\ll 1,$ (1)

where $b=\max(R_{1},R_{2},\Delta)$. The condition is necessary to employ the
static approach.

Figure 1: Sketch of geometry.

Consider the equation for potential $\phi$ under wave of $p$-polarization:

$\Delta\phi+k^{2}\phi=0.$ (2)

Then one can decompose potential $\phi$ into a series that includes only even
orders:

$\phi(x)=\phi^{(0)}(x)+k^{2}\phi^{(2)}(x)+\dots$ (3)

We get a chain of equations

$\displaystyle\Delta\phi^{(0)}(x)=0,$ (4)
$\displaystyle\Delta\phi^{(n)}(x)=-\phi^{(n-2)}(x),$ (5)

where $n\geqslant 1$.

Boundary conditions are

$\displaystyle\left.\varepsilon_{in}\frac{\partial\phi_{in}}{\partial{n}}\right|_{\gamma}=\left.\frac{\partial\phi_{out}}{\partial{n}}\right|_{\gamma},\left.\phi_{in}\right|_{\gamma}=\left.\phi_{out}\right|_{\gamma},$
(6)

where $\phi_{in}$ is the potential inside the cylinder, $\phi_{out}$ is the
potential outside, $\varepsilon_{in}$ is the dielectric permittivity, $\gamma$
is the boundary of the cylinder, $\vec{n}$ is the unit vector normal to the
border. So then, one has to solve the Laplace equation (4) in the zero-order.
Then, for the next orders we deal with Poisson equation (5). One can find all
next orders applying Green’s function for the Poisson equation derived in the
Appendix.

These conditions will apply to all orders, except of the zero-order, where one
has to include the incident field. Therefore, zero-order boundary conditions
are:

$\displaystyle\left.\varepsilon_{in}\frac{\partial\phi^{(0)}_{in}}{\partial\vec{n}}\right|_{\gamma}=\left.\frac{\partial\phi^{(0)}_{out}}{\partial\vec{n}}\right|_{\gamma}+\left.\frac{\partial\phi_{0}}{\partial
n}\right|_{\gamma},\quad\left.\phi^{(0)}_{in}\right|_{\gamma}=\left.\phi^{(0)}_{out}\right|_{\gamma}+\left.\phi_{0}\right|_{\gamma},$
(7)

where $\phi_{0}$ is the potential of incident field. In quasistatic
approximation the field is uniform, and then its potential is a linear
coordinate function

$\phi_{0}=-E_{0}y.$ (8)

##  III Conformal transformation

After a conformal mapping, a harmonic function remains harmonic ablowitz2003 ;
brown2009 . Then we look for a conformal transformation to build up an
analytical formula. To construct the map, we must analyze its level lines that
must walk along boundaries.

The appropriate map (bipolar coordinates) is described by the following
formula vorobev2010 :

$w(z)=\xi+i\eta=\ln\left(\frac{z+iC}{z-iC}\right),$ (9)

where $C$ is a real constant value. The transformation can be also rewritten
as

$e^{\xi}e^{i\eta}=\frac{x+iy+iC}{x+iy-iC}.$ (10)

Direct and inverse transformations are

$\displaystyle x=\frac{C\sin\eta}{\cosh\xi-\cos\eta},$ $\displaystyle
y=\frac{C\sinh\xi}{\cosh\xi-\cos\eta}.$
$\displaystyle\tanh\xi=\frac{2Cy}{x^{2}+y^{2}+C^{2}},$
$\displaystyle\tan\eta=\frac{2Cx}{x^{2}+y^{2}+C^{2}}.$ (11)

Excluding $\xi,\eta$ we can obtain parametric equations for level lines:

$\displaystyle
x^{2}+(y-C\coth\xi_{0})^{2}=\frac{C^{2}}{\sinh^{2}\xi_{0}},\quad\xi_{0}\neq
0,$
$\displaystyle(x+C\cot\eta_{0})^{2}+y^{2}=\frac{C^{2}}{\sin^{2}\eta_{0}},\quad\eta_{0}\neq
0.$ (12)

As one can see the first family of level lines are circles with radius
$|C/\sinh\xi_{0}|$ shifted to $\pm C\coth\xi_{0}$ in top and bottom semi-
planes. Similarly the second family consists of circles with radius
$|C/\sin\eta_{0}|$ shifted to $\pm C\cot\eta_{0}$ in left and right, as shown
in Fig. 2. Constant $C$ defines the coordinate origin and represents the
equipotential values for the borders. It also allows working with asymmetric
case, if we use top and bottom dashed circles of different radius. These
formulas define the origin of coordinate system and determine the boundary
values in new coordinates, when circle are mapped to constant $\xi$, with
$\eta\in(-\pi,\pi)$. The level lines are circles or their segments. The
transformation properly maps the boundaries of top and bottom circles into
equipotential lines.

Laplace’s operator in the $\eta$, $\xi$ coordinates is

$\Delta_{\xi,\eta}=\frac{(\cosh\xi-\cos\eta)^{2}}{C^{2}}\left(\frac{\partial^{2}}{\partial\xi^{2}}+\frac{\partial^{2}}{\partial\eta^{2}}\right).$
(13)

For circles $R_{1},R_{2}$ and gap dimension $\Delta$, one obtain the following
equations

$\displaystyle\frac{C}{\tanh\xi_{1}}=R_{1}+a_{1},\quad\frac{C}{\sinh\xi_{1}}=R_{1},$
$\displaystyle\frac{C}{\tanh\xi_{2}}=-R_{2}-a_{2},\quad\frac{C}{\sinh\xi_{2}}=-R_{2},$
$\displaystyle a_{1}+a_{2}=\Delta,$ (14)

where $\xi_{1},\xi_{2}$ are equipotential values for the boundaries. One can
see that Eq. (14) is symmetric with respect to transformation $y\to-y$:
$R_{1}\leftrightarrow-R_{2},a_{1}\leftrightarrow-a_{2},\Delta\to-\Delta$.

Solving equations we get

$\displaystyle\xi_{1}=\ln\left(1+\frac{a_{1}}{R_{1}}+\sqrt{\left(1+\frac{a_{1}}{R_{1}}\right)^{2}-1}\right),$
$\displaystyle\xi_{2}=\ln\left(1+\frac{a_{2}}{R_{2}}-\sqrt{\left(1+\frac{a_{2}}{R_{2}}\right)^{2}-1}\right),$
$\displaystyle a_{1}=\frac{\Delta(\Delta+2R_{2})}{2(\Delta+R_{1}+R_{2})},\quad
a_{2}=\frac{\Delta(\Delta+2R_{1})}{2(\Delta+R_{1}+R_{2})},$ $\displaystyle
C=\frac{\sqrt{\Delta(\Delta+2R_{1})(\Delta+2R_{2})(\Delta+2R_{1}+2R_{2})}}{2(\Delta+R_{1}+R_{2})}.$
(15)

In limiting cases of a narrow or wide slit the following asymptotics are valid

$C\approx\begin{cases}\sqrt{\frac{2\Delta
R_{1}R_{2}}{R_{1}+R_{2}}}+O(\Delta^{3/2}),&\Delta\ll R_{1,2},\\\
\frac{1}{2}(\Delta+R_{1}+R_{2})+O(\Delta^{-1}),&\Delta\gg R_{1,2}.\end{cases}$
(16)

The constant $C$ is plotted in Fig. 3 as a function of slit dimension
$\Delta$. We see the characteristic behaviour (16): at small $\Delta\ll
R_{1,2}$ the square root dependence takes place $C\sim\Delta^{1/2}$; at large
$\Delta\gg R_{1,2}$ it becomes linear $C\sim\Delta$.

Figure 2: Level lines (11) of the transformation in $(x,y)$ plane at $C=1$:
$\xi=0.5,1,1.5,2,\eta\in(-\pi,\pi)$ (solid thick lines);
$\xi=-0.5,-1,-1.5,-2,\eta\in(-\pi,\pi)$ (solid thin);
$\eta=0.5,1,1.5,2,\xi\in(-\pi,\pi)$ (dashed thick);
$\eta=-0.5,-1,-1.5,-2,\xi\in(-\pi,\pi)$ (dashed thin). Figure 3: Constant $C$
as a function of slit dimension $\Delta$ at fixed $R_{1}=0.1~{}\mu$m and
different $R_{2}=1,$ (dashed line), 0.5 (dot-dashed), 0.25 (thin), 0.125
(thick)

For two equal cylinders
$R_{1}=R_{2}=R,\varepsilon_{1}=\varepsilon_{2},a_{1}=a_{2}=a$ Eq. (15) yields
$a=\Delta/2,\xi_{1}=-\xi_{2},C\approx\sqrt{aR}$ ($a\ll R$). The solution turns
into symmetric problem studied by Vorob’ov vorobev2010 . At $\Delta,C\to 0$,
in the limit of ”kissing cylinders”, Eq. (9) reduces to an inversion
transformation introduced by Lei et al lei2010 :

$w(z)\approx\frac{2iC}{z}+O(C^{3}).$

Let us return to the general asymmetric case.

Boundary conditions (7) in bipolar coordinates take into account two borders:

$\displaystyle\left.\phi^{(0)}_{in}(x)\right|_{\xi_{1}}=\left.\phi^{(0)}_{out}(x)\right|_{\xi_{1}}+\left.\phi_{0}(x)\right|_{\xi_{1}},$
$\displaystyle\left.\phi^{(0)}_{in}(x)\right|_{\xi_{2}}=\left.\phi^{(0)}_{out}(x)\right|_{\xi_{2}}+\left.\phi_{0}(x)\right|_{\xi_{2}},$
$\displaystyle\varepsilon_{1}\left.\frac{\partial\phi^{(0)}_{in}(x)}{\partial\xi}\right|_{\xi_{1}}=\left.\frac{\partial\phi^{(0)}_{out}(x)}{\partial\xi}\right|_{\xi_{1}}+\left.\frac{\partial\phi_{0}(x)}{\partial\xi}\right|_{\xi_{1}},$
$\displaystyle\varepsilon_{2}\left.\frac{\partial\phi^{(0)}_{in}(x)}{\partial\xi}\right|_{\xi_{2}}=\left.\frac{\partial\phi^{(0)}_{out}(x)}{\partial\xi}\right|_{\xi_{2}}+\left.\frac{\partial\phi_{0}(x)}{\partial\xi}\right|_{\xi_{2}}.$
(17)

##  IV Quasistatic approximation

It is easy to find eigen functions of the Laplace operator in bipolar
coordinates

$\psi_{m}(\xi,\eta)=e^{\pm m\xi}\cos(m\eta).$ (18)

Then, we can determine the zero order as a series of the eigen functions:

$\phi^{(0)}=\begin{cases}\sum_{m=0}^{+\infty}A^{(0)}_{m}e^{-m\xi}\cos(m\eta),&\xi>\xi_{1},\\\
\sum_{m=0}^{+\infty}\left(C^{(0)}_{m}e^{m\xi}+D^{(0)}_{m}e^{-m\xi}\right)\cos(m\eta),&\xi_{2}<\xi<\xi_{1},\\\
\sum_{m=0}^{+\infty}B^{(0)}_{m}e^{m\xi}\cos(m\eta),&\xi<\xi_{2}.\end{cases}$
(19)

Area inside the first/second cylinder corresponds to
$\xi>\xi_{1}$/$\xi<\xi_{2}$, respectively. The interval $\xi_{2}<\xi<\xi_{1}$
is an external space. After solving boundary conditions equations at borders
$\xi=\xi_{1},\xi_{2}$ we obtain

$\displaystyle
A^{(0)}_{m}=\frac{-4E_{0}Ce^{2m\xi_{1}}\left(e^{2m\xi_{2}}(\varepsilon_{2}-1)+(\varepsilon_{2}+1)\right)}{J_{m}},$
$\displaystyle
B^{(0)}_{m}=\frac{4E_{0}C\left(e^{2m\xi_{1}}(\varepsilon_{1}+1)+(\varepsilon_{1}-1)\right)}{J_{m}},$
$\displaystyle
C^{(0)}_{m}=\frac{2E_{0}C(\varepsilon_{1}-1)\left(e^{2m\xi_{2}}(\varepsilon_{2}-1)+(\varepsilon_{2}+1)\right)}{J_{m}},$
$\displaystyle
D^{(0)}_{m}=\frac{-2E_{0}Ce^{2m\xi_{2}}(\varepsilon_{2}-1)\left(e^{2m\xi_{1}}(\varepsilon_{1}+1)+(\varepsilon_{1}-1)\right)}{J_{m}}.$
$\displaystyle
J_{m}=e^{2m\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2m\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1),\quad
m\neq 0.$ (20)

These formulas are appropriate to find the intensity distribution in the
cylinders vicinity. Returning to $(x,y)$ coordinates, we take the gradient to
determine the electric field. We keep the first 20 terms from the infinite
series as a trade-off between computational complexity and precision. We take
the amplitude of incident wave equal to 1. Fig. 4 and 5 present the intensity
distribution in arbitrary units. In both figures (a) and (c) are analytical
solutions, (b) and (d) are obtained by COMSOL Multiphysics. The level lines
demonstrate a good agreement of obtained formulas with the numerical
calculation.

We see the expected enhancement of electric field in the point between
circles. Local field enhancement by two cylinders due to their interaction and
strong optical coupling into a dimer are studied earlier in details
dmitriev2019 ; bulgakov2019 . For asymmetric situations the point of maximum
is shifted with respect to the geometric center.

Figure 4: Distribution of the intensity $|E|^{2}$ at $R_{1}=R_{2}=20$ nm,
$\varepsilon_{2}=2.25$, $\Delta=5\,$ nm, $\lambda=5\ \mu$m:
$\varepsilon_{1}=2.25$ for (a), (b) and $\varepsilon_{1}=2$ for (c), (d).

Figure 5: Distribution of the intensity $|E|^{2}$ at $R_{2}=20$ nm,
$\varepsilon_{1}=\varepsilon_{2}=2.25$, $\Delta=5\,$ nm, $\lambda=5\ \mu$m:
$R_{1}=40$ nm for (a), (b) and $R_{1}=60$ nm for (c), (d).

##  V Conclusions

We obtain the quasistatic approximation of $p$-wave scattered on two parallel
cylinders. For this purpose a conformal transformation in the plane of cross-
section is built up that maps the borders into level lines of the new
orthogonal coordinate system, the bipolar coordinates. In the symmetric case,
where cylinders are of equal diameters and refractive indices, the solution is
shown to tend to known formulas. The field can be expressed as a decomposition
in terms of the Laplace operator eigenfunctions. The system of algebraic
equations for its coefficients has been solved and the near-field intensity
distribution has been found. This solution agrees with the numerical
calculations by COMSOL Multiphysics. Finally, we build up Green’s function,
making it possible to calculate the higher orders.

The formulas obtained are useful for further investigation of the scattering
process, in particular, for correct calculation of resonant configurations in
a nanoparticle structure. The quasistatic approximation is valid while the
diameters and gap between the cylinders remain much less than the radiation
wavelength.

## Acknowledgements

The work was supported by Russian Science Foundation, grant #22-22-00633.

## Appendix: Green’s function in bipolar coordinates

Similar to Eq. (19), we can write the Green function as a series of the eigen
functions. In the bipolar coordinates free space is mapped into the strip
$\xi_{2}<\xi<\xi_{1}$. There are three epressions of the function in different
domains:

$G(\xi,\xi^{\prime},\eta,\eta^{\prime})=\begin{cases}G_{1}(\xi,\xi^{\prime},\eta,\eta^{\prime}),\xi^{\prime}>\xi_{1}\\\
G_{2}(\xi,\xi^{\prime},\eta,\eta^{\prime}),\xi_{1}>\xi^{\prime}>\xi_{2}\\\
G_{3}(\xi,\xi^{\prime},\eta,\eta^{\prime}),\xi_{2}>\xi^{\prime}\end{cases}$
(A.1)

In the first domain

$G_{1}(\xi,\xi^{\prime},\eta,\eta^{\prime})=\begin{cases}\sum_{m=0}^{+\infty}A^{1}_{m}e^{-m\xi}\cos(m(\eta-\eta^{\prime})),\xi>\xi^{\prime},\\\
\sum_{m=0}^{+\infty}\left(F^{1}_{m}e^{m\xi}+L^{1}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi^{\prime}>\xi>\xi_{1},\\\
\sum_{m=0}^{+\infty}\left(C^{1}_{m}e^{m\xi}+D^{1}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi_{1}>\xi>\xi_{2},\\\
\sum_{m=0}^{+\infty}B^{1}_{m}e^{m\xi}\cos(m(\eta-\eta^{\prime})),\xi>\xi^{\prime},\end{cases}$
(A.2)

where

$\displaystyle\begin{gathered}A^{1}_{n}=-\frac{\alpha
e^{-n\xi^{\prime}}\left(e^{4n\xi_{1}}(\varepsilon_{1}-1)(\varepsilon_{2}+1)+e^{2n(\xi_{1}+\xi^{\prime})}(\varepsilon_{1}+1)(\varepsilon_{2}+1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)}+\\\
+\frac{\alpha
e^{-n\xi^{\prime}}\left(e^{2n(\xi_{2}+\xi^{\prime})}(\varepsilon_{1}-1)(\varepsilon_{2}-1)+e^{2n(\xi_{1}+\xi_{2})}(\varepsilon_{1}+1)(\varepsilon_{2}-1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},\end{gathered}$
$\displaystyle
B^{1}_{n}=-\frac{2\alpha\varepsilon_{1}e^{2n\xi_{1}-n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle
C^{1}_{n}=-\frac{\alpha\varepsilon_{1}(\varepsilon_{2}+1)e^{2n\xi_{1}-n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle
D^{1}_{n}=-\frac{\alpha\varepsilon_{1}(\varepsilon_{2}-1)e^{2n\xi_{1}+2n\xi_{2}-n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle L^{1}_{n}=-\frac{\alpha
e^{2n\xi_{1}-n\xi^{\prime}}\left(e^{2n\xi_{1}}(\varepsilon_{1}-1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}+1)(\varepsilon_{2}-1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},$
$\displaystyle F^{1}_{n}=-\frac{1}{2}\alpha
e^{-n\xi^{\prime}},\quad\alpha=\frac{C^{2}}{\left(\cosh\xi^{\prime}-\cos\eta^{\prime}\right)^{2}}.$

The second function is

$G_{2}(\xi,\xi^{\prime},\eta,\eta^{\prime})=\begin{cases}\sum_{m=0}^{+\infty}A^{2}_{m}e^{-m\xi}\cos(m(\eta-\eta^{\prime})),\xi>\xi_{1},\\\
\sum_{m=0}^{+\infty}\left(C^{2}_{m}e^{m\xi}+D^{2}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi_{1}>\xi>\xi^{\prime},\\\
\sum_{m=0}^{+\infty}\left(F^{2}_{m}e^{m\xi}+L^{2}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi^{\prime}>\xi>\xi_{2},\\\
\sum_{m=0}^{+\infty}B^{2}_{m}e^{m\xi}\cos(m(\eta-\eta^{\prime})),\xi<\xi_{2},\end{cases}$
(A.4)

where

$\displaystyle A^{2}_{n}=-\frac{\alpha
e^{2n\xi_{1}-n\xi^{\prime}}\left(-e^{2n\xi_{2}}(\varepsilon_{2}-1)+e^{2n\xi^{\prime}}(\varepsilon_{2}+1)\right)}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle B^{2}_{n}=\frac{\alpha
e^{-n\xi_{3}}\left(e^{2n\xi^{\prime}}(\varepsilon_{1}-1)-e^{2n\xi_{1}}(\varepsilon_{1}+1)\right)}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle C^{2}_{n}=\frac{\alpha
e^{-n\xi_{3}}(\varepsilon_{1}-1)\left(e^{2n\xi^{\prime}}(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{2}-1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},$
$\displaystyle D^{2}_{n}=-\frac{\alpha
e^{2n\xi_{1}-n\xi_{3}}(\varepsilon_{1}+1)\left(e^{2n\xi^{\prime}}(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{2}-1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},$
$\displaystyle F^{2}_{n}=-\frac{\alpha
e^{-n\xi_{3}}(\varepsilon_{2}+1)\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)-e^{2n\xi^{\prime}}(\varepsilon_{1}-1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},$
$\displaystyle L^{2}_{n}=-\frac{\alpha
e^{2n\xi_{2}-n\xi_{3}}(\varepsilon_{2}-1)\left(e^{2n\xi^{\prime}}(\varepsilon_{1}-1)-e^{2n\xi_{1}}(\varepsilon_{1}+1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)}.$

The third value is

$G_{3}(\xi,\xi^{\prime},\eta,\eta^{\prime})=\begin{cases}\sum_{m=0}^{+\infty}A^{3}_{m}e^{-m\xi}\cos(m(\eta-\eta^{\prime})),\xi>\xi_{1},\\\
\sum_{m=0}^{+\infty}\left(C^{3}_{m}e^{m\xi}+D^{3}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi_{1}>\xi>\xi_{2},\\\
\sum_{m=0}^{+\infty}\left(F^{3}_{m}e^{m\xi}+L^{3}_{m}e^{-m\xi}\right)\cos(m(\eta-\eta^{\prime})),\xi_{2}>\xi>\xi^{\prime},\\\
\sum_{m=0}^{+\infty}B^{3}_{m}e^{m\xi}\cos(m(\eta-\eta^{\prime})),\xi^{\prime}>\xi,\end{cases},$
(A.5)

where

$\displaystyle
A^{3}_{n}=-\frac{2\alpha\varepsilon_{2}e^{2n\xi_{1}+n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle\begin{gathered}B^{3}_{n}=-\frac{\alpha
e^{-n(2\xi_{2}+\xi^{\prime})}\left(e^{4n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)+e^{2n(\xi_{2}+\xi^{\prime})}(\varepsilon_{1}-1)(\varepsilon_{2}+1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)}+\\\
+\frac{\alpha
e^{-n(2\xi_{2}+\xi^{\prime})}\left(e^{2n(\xi_{1}+\xi^{\prime})}(\varepsilon_{1}+1)(\varepsilon_{2}-1)+e^{2n(\xi_{1}+\xi_{2})}(\varepsilon_{1}+1)(\varepsilon_{2}+1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},\end{gathered}$
$\displaystyle
C^{3}_{n}=\frac{\alpha\varepsilon_{2}(\varepsilon_{1}-1)e^{n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle
D^{3}_{n}=-\frac{\alpha\varepsilon_{2}(\varepsilon_{1}+1)e^{2n\xi_{1}+n\xi^{\prime}}}{e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)},$
$\displaystyle F^{3}_{n}=-\frac{\alpha
e^{-2n\xi_{2}+n\xi^{\prime}}\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}-1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}+1)\right)}{2\left(e^{2n\xi_{1}}(\varepsilon_{1}+1)(\varepsilon_{2}+1)-e^{2n\xi_{2}}(\varepsilon_{1}-1)(\varepsilon_{2}-1)\right)},L^{3}_{n}=-\frac{1}{2}\alpha
e^{n\xi^{\prime}}.$

## References

  * (1) C. Girard, “Near fields in nanostructures,” _Reports on progress in physics_ , vol. 68, no. 8, p. 1883, 2005.
  * (2) V. M. Shalaev and S. Kawata, _Nanophotonics with surface plasmons_. Elsevier, 2006.
  * (3) M. I. Stockman, “Nanoplasmonics: The physics behind the applications,” _Phys. Today_ , vol. 64, no. 2, pp. 39–44, 2011.
  * (4) L. Novotny and B. Hecht, _Principles of nano-optics_. Cambridge university press, 2012.
  * (5) A. Devilez, X. Zambrana-Puyalto, B. Stout, and N. Bonod, “Mimicking localized surface plasmons with dielectric particles,” _Physical Review B_ , vol. 92, no. 24, p. 241412, 2015.
  * (6) A. Krasnok, M. Caldarola, N. Bonod, and A. Alú, “Spectroscopy and biosensing with optically resonant dielectric nanostructures,” _Advanced optical materials_ , vol. 6, no. 5, p. 1701094, 2018.
  * (7) F. Enrichi, A. Quandt, and G. C. Righini, “Plasmonic enhanced solar cells: Summary of possible strategies and recent results,” _Renewable and Sustainable Energy Reviews_ , vol. 82, pp. 2433–2439, 2018.
  * (8) J. Homola, “Surface plasmon resonance sensors for detection of chemical and biological species,” _Chemical Reviews_ , vol. 108, no. 2, pp. 462–493, 2008\.
  * (9) A. Shalabney and I. Abdulhalim, “Sensitivity-enhancement methods for surface plasmon sensors,” _Laser & Photonics Reviews_, vol. 5, no. 4, pp. 571–606, 2011.
  * (10) M. Piliarik, H. Šípová, P. Kvasnička, N. Galler, J. R. Krenn, and J. Homola, “High-resolution biosensor based on localized surface plasmons,” _Opt. Express_ , vol. 20, no. 1, pp. 672–680, 2012.
  * (11) Y. Xu, P. Bai, X. Zhou, Y. Akimov, C. E. Png, L.-K. Ang, W. Knoll, and L. Wu, “Optical refractive index sensors with plasmonic and photonic structures: promising and inconvenient truth,” _Advanced Optical Materials_ , vol. 7, no. 9, p. 1801433, 2019.
  * (12) A. Bialiayeu, A. Bottomley, D. Prezgot, A. Ianoul, and J. Albert, “Plasmon-enhanced refractometry using silver nanowire coatings on tilted fibre bragg gratings,” _Nanotechnology_ , vol. 23, no. 44, p. 444012, 2012\.
  * (13) P. T. Arasu, A. S. M. Noor, A. A. Shabaneh, M. H. Yaacob, H. N. Lim, and M. A. Mahdi, “Fiber bragg grating assisted surface plasmon resonance sensor with graphene oxide sensing layer,” _Opt. Commun._ , vol. 380, pp. 260–266, 2016\.
  * (14) K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” _IEEE Transactions on antennas and propagation_ , vol. 14, no. 3, pp. 302–307, 1966.
  * (15) E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” _The Astrophysical Journal_ , vol. 186, pp. 705–714, 1973.
  * (16) B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for periodic targets: theory and tests,” _J. Optical Society of America A_ , vol. 25, no. 11, pp. 2693–2703, 2008.
  * (17) C. A. Brebbia and S. Walker, _Boundary element techniques in engineering_. Elsevier, 2016.
  * (18) O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, _The finite element method: its basis and fundamentals_. Elsevier, 2005.
  * (19) N. N. Voitovich, B. Z. Katsenelenbaum, and A. N. Sivov, “The generalized natural-oscillation method in diffraction theory,” _Sov. Phys. Usp._ , vol. 19, no. 9, pp. 337–352, 1976.
  * (20) P. Vorobev, “Electric field enhancement between two parallel cylinders due to plasmonic resonance,” _Journal of Experimental and Theoretical Physics_ , vol. 110, no. 2, pp. 193–198, 2010.
  * (21) D. Y. Lei, A. Aubry, S. A. Maier, and J. B. Pendry, “Broadband nano-focusing of light using kissing nanowires,” _New Journal of Physics_ , vol. 12, no. 9, p. 093030, 2010.
  * (22) A. Bereza, A. Nemykin, S. Perminov, L. Frumin, and D. Shapiro, “Light scattering by dielectric bodies in the born approximation,” _Physical Review A_ , vol. 95, no. 6, p. 063839, 2017.
  * (23) A. Bereza, L. Frumin, A. Nemykin, S. Perminov, and D. Shapiro, “Perturbation series for the scattering of electromagnetic waves by parallel cylinders,” _EPL (Europhysics Letters)_ , vol. 127, no. 2, p. 20002, 2019.
  * (24) T. A. van der Sijs, O. El Gawhary, and H. P. Urbach, “Electromagnetic scattering beyond the weak regime: Solving the problem of divergent born perturbation series by padé approximants,” _Phys. Rev. Research_ , vol. 2, p. 013308, Mar 2020.
  * (25) J. A. Rebouças and P. A. Brandão, “Scattering of light by a parity-time-symmetric dipole beyond the first born approximation,” _Phys. Rev. A_ , vol. 104, p. 063514, Dec 2021.
  * (26) N. A. Ustimenko, D. F. Kornovan, K. V. Baryshnikova, A. B. Evlyukhin, and M. I. Petrov, “Multipole born series approach to light scattering by mie-resonant nanoparticle structures,” _Journal of Optics_ , vol. 24, no. 3, p. 035603, 2022.
  * (27) S. Belan and S. Vergeles, “Plasmon mode propagation in array of closely spaced metallic cylinders,” _Optical Materials Express_ , vol. 5, no. 1, pp. 130–141, 2015.
  * (28) S.-C. Lee, “Scattering at oblique incidence by multiple cylinders in front of a surface,” _Journal of Quantitative Spectroscopy and Radiative Transfer_ , vol. 182, pp. 119–127, 2016.
  * (29) M. J. Ablowitz, A. S. Fokas, and A. S. Fokas, _Complex variables: introduction and applications_. Cambridge University Press, 2003.
  * (30) J. W. Brown and R. V. Churchill, _Complex variables and applications_. McGraw-Hill, 2009.
  * (31) A. A. Dmitriev and M. V. Rybin, “Combining isolated scatterers into a dimer by strong optical coupling,” _Physical Review A_ , vol. 99, no. 6, p. 063837, 2019.
  * (32) E. N. Bulgakov, K. N. Pichugin, and A. F. Sadreev, “Evolution of the resonances of two parallel dielectric cylinders with distance between them,” _Physical Review A_ , vol. 100, no. 4, p. 043806, 2019.


* Federrath et al. (2011a) Federrath C., Banerjee R., Seifried D., Clark P. C., Klessen R. S., 2011a, in J. Alves, B. G. Elmegreen, J. M. Girart, & V. Trimble ed., IAU Symposium Vol. 270, Computational Star Formation. pp 425–428 (arXiv:1007.2504), doi:10.1017/S1743921311000755
  * Federrath et al. (2011b) Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011b, ApJ, 731, 62
  * Federrath et al. (2014) Federrath C., Schrön M., Banerjee R., Klessen R. S., 2014, ApJ, 790, 128
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy, 5, 365
  * Federrath et al. (2022) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2022, TG: Turbulence Generator, Astrophysics Source Code Library, record ascl:2204.001 (ascl:2204.001)
  * Finn et al. (2019) Finn M. K., Johnson K. E., Brogan C. L., Wilson C. D., Indebetouw R., Harris W. E., Kamenetzky J., Bemis A., 2019, ApJ, 874, 120
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Fukushima & Yajima (2021) Fukushima H., Yajima H., 2021, MNRAS, 506, 5512
  * Gallegos-Garcia et al. (2020) Gallegos-Garcia M., Burkhart B., Rosen A. L., Naiman J. P., Ramirez-Ruiz E., 2020, ApJ, 899, L30
  * Gavagnin et al. (2017) Gavagnin E., Bleuler A., Rosdahl J., Teyssier R., 2017, MNRAS, 472, 4155
  * Geen et al. (2016) Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2016, MNRAS, 463, 3129
  * Gentry et al. (2019) Gentry E. S., Krumholz M. R., Madau P., Lupi A., 2019, MNRAS, 483, 3647
  * Grasha et al. (2019) Grasha K., et al., 2019, MNRAS, 483, 4707
  * Gritschneder et al. (2009) Gritschneder M., Naab T., Walch S., Burkert A., Heitsch F., 2009, Astrophysical Journal, 694
  * Grudić et al. (2018) Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, MNRAS, 475, 3511
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * He et al. (2020) He H., Wilson C. D., Sliwa K., Iono D., Saito T., 2020, MNRAS, 496, 5243
  * He et al. (2022) He H., Wilson C., Brunetti N., Finn M., Bemis A., Johnson K., 2022, ApJ, 928, 57
  * Heyer et al. (2016) Heyer M., Gutermuth R., Urquhart J. S., Csengeri T., Wienen M., Leurini S., Menten K., Wyrowski F., 2016, A&A, 588, A29
  * Hollyhead et al. (2015) Hollyhead K., Bastian N., Adamo A., Silva-Villa E., Dale J., Ryon J. E., Gazak Z., 2015, MNRAS, 449, 1106
  * Hopkins (2018) Hopkins A. M., 2018, preprint,  (arXiv:1807.09949)
  * Hopkins et al. (2010) Hopkins P. F., Murray N., Quataert E., Thompson T. A., 2010, MNRAS, 401, L19
  * Hopkins et al. (2011) Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3522
  * Hu et al. (2022) Hu Z., Krumholz M. R., Pokhrel R., Gutermuth R. A., 2022, MNRAS, 511, 1431
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Jacquet & Krumholz (2011) Jacquet E., Krumholz M. R., 2011, ApJ, 730, 116
  * Jiang (2021) Jiang Y.-F., 2021, ApJS, 253, 49
  * Jiang et al. (2013) Jiang Y.-F., Davis S. W., Stone J. M., 2013, ApJ, 763, 102
  * Kim & Ostriker (2017) Kim C.-G., Ostriker E. C., 2017, ApJ, 846, 133
  * Kim et al. (2017) Kim J.-G., Kim W.-T., Ostriker E. C., Skinner M. A., 2017, ApJ, 851, 93
  * Kim et al. (2018) Kim J.-G., Kim W.-T., Ostriker E. C., 2018, ApJ, 859, 68
  * Krumholz (2018) Krumholz M. R., 2018, MNRAS, 480, 3468
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352
  * Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304
  * Krumholz & Thompson (2012) Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155
  * Krumholz & Thompson (2013) Krumholz M. R., Thompson T. A., 2013, Monthly Notices of the Royal Astronomical Society, 434, 2329
  * Krumholz et al. (2006) Krumholz M. R., Matzner C. D., McKee C. F., 2006, ApJ, 653, 361
  * Krumholz et al. (2007) Krumholz M. R., Klein R. I., McKee C. F., Bolstad J., 2007, ApJ, 667, 626
  * Krumholz et al. (2009) Krumholz M. R., Klein R. I., McKee C. F., Offner S. S. R., Cunningham A. J., 2009, Science, 323, 754
  * Krumholz et al. (2012a) Krumholz M. R., Dekel A., McKee C. F., 2012a, ApJ, 745, 69
  * Krumholz et al. (2012b) Krumholz M. R., Klein R. I., McKee C. F., 2012b, ApJ, 754, 71
  * Krumholz et al. (2014) Krumholz M. R., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 243 (arXiv:1401.2473), doi:10.2458/azu_uapress_9780816531240-ch011
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227
  * Kuiper et al. (2012) Kuiper R., Klahr H., Beuther H., Henning T., 2012, A&A, 537, A122
  * Lancaster et al. (2021a) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89
  * Lancaster et al. (2021b) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 922, L3
  * Lane et al. (2022) Lane H. B., Grudić M. Y., Guszejnov D., Offner S. S. R., Faucher-Giguère C.-A., Rosen A. L., 2022, MNRAS, 510, 4767
  * Leroy et al. (2017) Leroy A. K., et al., 2017, ApJ, 846, 71
  * Leroy et al. (2018) Leroy A. K., et al., 2018, The Astrophysical Journal, 869, 126
  * Levy et al. (2021) Levy R. C., et al., 2021, ApJ, 912, 4
  * Lopez et al. (2011) Lopez L. A., Krumholz M. R., Bolatto A. D., Prochaska J. X., Ramirez-Ruiz E., 2011, ApJ, 731, 91
  * Lopez et al. (2014) Lopez L. A., Krumholz M. R., Bolatto A. D., Prochaska J. X., Ramirez-Ruiz E., Castro D., 2014, ApJ, 795, 121
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * MacNeice et al. (2000) MacNeice P., Olson K. M., Mobarry C., de Fainchtein R., Packer C., 2000, Computer Physics Communications, 126, 330
  * Markevitch & Vikhlinin (2007) Markevitch M., Vikhlinin A., 2007, Phys. Rep., 443, 1
  * Marshak (1958) Marshak R. E., 1958, Physics of Fluids, 1, 24
  * McCrady et al. (2005) McCrady N., Graham J. R., Vacca W. D., 2005, ApJ, 621, 278
  * Meidt et al. (2018) Meidt S. E., et al., 2018, ApJ, 854, 100
  * Menon et al. (2020) Menon S. H., Federrath C., Kuiper R., 2020, MNRAS, 493, 4643
  * Menon et al. (2021) Menon S. H., Federrath C., Klaassen P., Kuiper R., Reiter M., 2021, MNRAS, 500, 1721
  * Menon et al. (2022) Menon S. H., Federrath C., Krumholz M. R., Kuiper R., Wibking B. D., Jung M., 2022, MNRAS, 512, 401
  * Mihalas & Auer (2001) Mihalas D., Auer L., 2001, Journal of Quantitative Spectroscopy and Radiative Transfer, 71, 61
  * Mihalas & Klein (1982) Mihalas D., Klein R. I., 1982, Journal of Computational Physics, 46, 97
  * Murray et al. (2010a) Murray N., Quataert E., Thompson T. A., 2010a, The Astrophysical Journal, 709, 191
  * Murray et al. (2010b) Murray N., Quataert E., Thompson T. A., 2010b, ApJ, 709, 191
  * Nakamura & Li (2007) Nakamura F., Li Z.-Y., 2007, ApJ, 662, 395
  * Offner & Liu (2018) Offner S. S. R., Liu Y., 2018, Nature Astronomy, 2, 896
  * Olivier et al. (2021) Olivier G. M., Lopez L. A., Rosen A. L., Nayak O., Reiter M., Krumholz M. R., Bolatto A. D., 2021, ApJ, 908, 68
  * Ostriker et al. (2010) Ostriker E. C., McKee C. F., Leroy A. K., 2010, ApJ, 721, 975
  * Padoan et al. (2014) Padoan P., Federrath C., Chabrier G., Evans II N. J., Johnstone D., Jørgensen J. K., McKee C. F., Nordlund Å., 2014, in , Protostars and Planets VI. University of Arizona Press (arXiv:1312.5365), doi:10.2458/azu_uapress_9780816531240-ch004
  * Padoan et al. (2016) Padoan P., Pan L., Haugbølle T., Nordlund Å., 2016, ApJ, 822, 11
  * Park et al. (2012) Park H., Knoll D. A., Rauenzahn R. M., Wollaber A. B., Densmore J. D., 2012, Transport Theory and Statistical Physics, 41, 284
  * Pokhrel et al. (2020) Pokhrel R., et al., 2020, ApJ, 896, 60
  * Portegies Zwart et al. (2010) Portegies Zwart S., McMillan S., Gieles M., 2010, Annual Review of Astronomy and Astrophysics, 48, 431
  * Raskutti et al. (2016) Raskutti S., Ostriker E. C., Skinner M. A., 2016, ApJ, 829, 130
  * Reissl et al. (2018) Reissl S., Klessen R. S., Low M. M. M., Pellegrini E. W., 2018, Astronomy and Astrophysics, 611
  * Ricker (2008) Ricker P. M., 2008, ApJS, 176, 293
  * Rico-Villas et al. (2020) Rico-Villas F., Martín-Pintado J., González-Alfonso E., Martín S., Rivilla V. M., 2020, MNRAS, 491, 4573
  * Rogers & Pittard (2013) Rogers H., Pittard J. M., 2013, MNRAS, 431, 1337
  * Rosdahl & Teyssier (2015) Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380
  * Rosen et al. (2014) Rosen A. L., Lopez L. A., Krumholz M. R., Ramirez-Ruiz E., 2014, MNRAS, 442, 2701
  * Rosen et al. (2016) Rosen A. L., Krumholz M. R., McKee C. F., Klein R. I., 2016, MNRAS, 463, 2553
  * Saad (2003) Saad Y., 2003, Iterative Methods for Sparse Linear Systems: Second Edition. Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, https://books.google.com.au/books?id=qtzmkzzqFmcC
  * Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  * Skinner & Ostriker (2013) Skinner M. A., Ostriker E. C., 2013, Astrophysical Journal, Supplement Series, 206
  * Skinner & Ostriker (2015) Skinner M. A., Ostriker E. C., 2015, ApJ, 809, 187
  * Smith et al. (2020) Smith L. J., Bajaj V., Ryon J., Sabbi E., 2020, ApJ, 896, 84
  * Thompson & Krumholz (2016) Thompson T. A., Krumholz M. R., 2016, Monthly Notices of the Royal Astronomical Society, 455, 334
  * Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167
  * Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman John H. I., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179
  * Tsang & Milosavljević (2015) Tsang B. T. H., Milosavljević M., 2015, MNRAS, 453, 1108
  * Tsang & Milosavljević (2018) Tsang B. T. H., Milosavljević M., 2018, MNRAS, 478, 4142
  * Turk et al. (2010) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2010, The Astrophysical Journal Supplement Series, 192, 9
  * Turner et al. (2015) Turner J. L., Beck S. C., Benford D. J., Consiglio S. M., Ho P. T. P., Kovács A., Meier D. S., Zhao J. H., 2015, Nature, 519, 331
  * Turner et al. (2017) Turner J. L., Consiglio S. M., Beck S. C., Goss W. M., Ho P. T. P., Meier D. S., Silich S., Zhao J.-H., 2017, ApJ, 846, 73
  * Utomo et al. (2018) Utomo D., et al., 2018, ApJ, 861, L18
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Waagan et al. (2011) Waagan K., Federrath C., Klingenberg C., 2011, Journal of Computational Physics, 230, 3331
  * Wang et al. (2010) Wang P., Li Z., Abel T., Nakamura F., 2010, ApJ, 709, 27
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Wibking et al. (2018) Wibking B. D., Thompson T. A., Krumholz M. R., 2018, MNRAS, 477, 4665
  * Zhang (2018) Zhang D., 2018, Galaxies, 6, 114
  * Zhang & Davis (2017) Zhang D., Davis S. W., 2017, ApJ, 839, 54
  * Zhang et al. (2018) Zhang D., Davis S. W., Jiang Y.-F., Stone J. M., 2018, ApJ, 854, 110

##  Appendix A Numerical convergence studies

In this section, we test the dependencies of our results on certain numerical
choices that have been made in our simulations. Since the fiducial simulations
are all sub-Eddington, the values of $\epsilon_{*}$ attained are virtually
identical, irrespective of these numerical choices. Therefore, we found it
more informative to compare the Eddington ratio $f_{\mathrm{Edd}}$ with
different numerical choices.

###  A.1 Dependence on source size

Figure 20: Comparison of the obtained Eddington ratios for the run S5KSemA2F1
at $t=3\,t_{\mathrm{ff}}$, with different choices for the source size
($\sigma_{*}$). The choice adopted in our study is $\sigma_{*}=2.5\Delta x$,
where $\Delta x$ is the cell width.

Here we test the dependence of our results on the adopted choice for the size
of the source $\sigma_{*}$ in Equation (15). Our fiducial choice is
$\sigma_{*}=2.5\Delta x$, where $\Delta x$ is the size of a computational cell
on the highest level of AMR where the sink particle is located; for our
fiducial resolution, this corresponds to $\sigma_{*}=0.039\,R_{\rm cloud}$. We
check the dependence of our results on this choice by repeating run S5KSemA2F1
with alternate choices of $\sigma_{*}=5\Delta x=0.078\,R_{\rm cloud}$ and
$\sigma_{*}=10\Delta x=0.16\,R_{\rm cloud}$. We plot the Eddington ratios we
obtain at $t=3t_{\mathrm{ff}}$ for this test in Figure 20. We see that all the
curves are nearly identical at all radii substantially larger than
$\sigma_{*}$; at $r\lesssim\sigma_{*}$, the Eddington ratio $f_{\mathrm{Edd}}$
decreases because $j_{*}$ is smoothed in this region, and hence radiation
forces are smaller. Based on this test we conclude that the effects of
smoothing are minimal for our fiducial choice of $\sigma_{*}$, because Figure
20 demonstrates that smoothing has negligible effects outside the smoothed
volume, and the volume that is affected for $\sigma_{*}=2.5\Delta x$ is only
$\sim 10^{-5}$ of the total cloud volume.

###  A.2 Dependence on the micro-physical dust model

Figure 21: Eddington ratios at $t=3\,t_{\mathrm{ff}}$ in run S5KSemA2F1
repeated using different micro-physical properties. The opacities in the
fiducial run assume that dust grains are aggregate particles made of composite
materials with normal iron content in the silicates. The alternative models
shown use spherical (Sph) rather then aggregate particles, homogeneous (Hom)
rather than composite dust, and iron-poor silicates (IPS) rather than normal
iron content. See text for further details.

Since we find that realistic dust opacities (i.e., $\kappa_{\mathrm{Sem}}$)
are crucial to our results, it is important to check for any dependencies on
the particulars of the dust model. Semenov et al. (2003) present a range of
possible models, based on differing assumptions regarding the micro-physical
nature of dust grains – specifically their topology, shape and chemical
composition. Our default choices are grains with an “aggregate” shape (i.e.,
they are a cluster of small spherical sub-grains sticking together),
composite” composition (i.e., they incorporate a mixture of materials), and
“normal” (NRM) iron content in the dust silicates; alternative assumptions for
each of these properties lead to changes in the opacity (see fig. 1 in Semenov
et al., 2003).

To test the sensitivity of our results to these choices, we again repeat run
S5KSemA2F1 using some of Semenov et al.’s alternative dust models. We test the
following choices: “spherical” particles instead of a cluster of sub-grains,
“homogeneous” rather than heterogeneous composition, and “Iron-poor silicates”
(IPS) instead of normal (NRM) iron content; for each of these variations, we
leave the remaining properties unchanged, e.g., our spherical grain models
have composite composition and normal iron content. We plot the Eddington
ratios $f_{\mathrm{Edd}}$ obtained from these experiments in Figure 21. We see
that alternative dust models lead to mild changes in the Eddington ratio, but
the results are qualitatively quite similar (the same to within a factor of
$\lesssim 2$), and most importantly, all are sub-Eddington. This suggests that
our major results are robust against plausible changes in the assumed micro-
physical dust properties.

###  A.3 Dependence on resolution

Figure 22: Comparison of the obtained Eddington ratios for the fiducial run at
$t=3\,t_{\mathrm{ff}}$, with different grid resolutions. Solid lines indicate
the resolution used in exploring parameter space in our study. The AMR run had
an effective maximum resolution of $1024^{3}$ cells. We see that the runs are
reasonably converged with resolution.

We test for numerical convergence of our results by comparing runs with
different grid resolutions. We again repeat simulation S5KSemA2F1 with
uniform-grid resolutions of $64^{3}$ and $128^{3}$ to compare with our
fiducial choice of $256^{3}$, and we further compare to an AMR run with a
base-grid resolution of $128^{3}$, and 3 levels of refinement, thereby leading
to a maximum effective resolution of $1024^{3}$. For the AMR run we refine
based on the Jeans length, requiring that it be resolved by at least 8 cells
(Truelove et al., 1997; Federrath et al., 2011b). We compare the radial
profiles of Eddington ratios in these runs in Figure 22. We see that at higher
resolution, the radial location of maximum $f_{\mathrm{Edd}}$ shifts inwards
with increasing resolution. However, the maximum value of $f_{\mathrm{Edd}}$
shows changes of at most $\sim 20\%$, and over the bulk of the cloud,
$f_{\mathrm{Edd}}$ is fairly independent of resolution. Therefore, we conclude
that our qualitative results are insensitive to resolution, and thereby are
reasonably converged.

##  Appendix B Idealised spherically-symmetric model cloud

Figure 23: The reduced flux ($\zeta$; top), radiation temperature
($T_{\mathrm{r}}$; middle), and Rosseland opacity ($\kappa_{\mathrm{R}}$;
bottom) profiles obtained with the steady-state semi-analytic models described
in Section B. Colours indicate the four values of $\Sigma_{\mathrm{cloud}}$
probed in our study; solid, dashed and dash-dotted lines indicate values of
$\alpha$ in the density profile.

Here we fully describe the semi-analytical model that we use to obtain the
steady-state radial profiles of matter and radiation quantities in spherical
symmetry. Our approach follows the one described in section 4.3.6 of Skinner &
Ostriker (2013), modified to accommodate non-constant opacities. We model a
spherically-symmetric cloud of gas and stars with a total mass
$M=10^{6}\,\mathrm{M}_{\sun}$; the stars are located at the origin ($r=0$),
and we inject (IR) photons with a radial distribution that follows Equation
(15). Their luminosity $L_{*}=\langle L_{*}/M_{*}\rangle M_{*}$, where
$M_{*}=\epsilon_{*}M$ is the stellar mass, $\epsilon_{*}$ is the mass fraction
in stars, and $\langle L_{*}/M_{*}\rangle=1.7\times
10^{3}\,\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{g}^{-1}$. We set the size of
the source $\sigma_{*}=0.03\,R_{\mathrm{g}}$, where $R_{\mathrm{g}}$ is the
maximum radius of the cloud, but, as discussed in Appendix A.1, this choice
makes little difference at radii $r>\sigma_{*}$. We parameterise our model in
terms of the gas mass surface density $\Sigma=M_{\mathrm{g}}/(\pi
R_{\mathrm{g}}^{2})$, and an index $\alpha$ such that $\rho\propto
r^{-\alpha}$ describing the radial profile of gas density. For a given
$\Sigma$ and $\alpha$, it is straightforward to obtain the normalisation of
the density profile required to reach the desired mass $M$. We set
$\rho(r)=10^{-10}\rho_{0}$ for $r>R_{\mathrm{g}}$, where
$\rho_{0}=M_{\mathrm{g}}/(4/3\pi R_{\mathrm{g}}^{3})$, to emulate an
optically-thin background outside the cloud.

Given the parameters, the radial radiation flux $F$ can then be obtained by
solving the steady-state equation of $E_{r}$ (Equation 5),

$\nabla\cdot F=j_{*}(r).$ (29)

This gives

$F_{*}(r)=\frac{L_{*}}{4\pi
r^{2}}\left[\operatorname{erf}\left(\frac{r}{\sqrt{2}\sigma_{*}}\right)-\frac{2r}{\sqrt{2\pi\sigma_{*}^{2}}}\exp\left(-\frac{r^{2}}{2\sigma_{*}^{2}}\right)\right].$
(30)

The equation for the radiation momentum (Equation 6) in steady state
satisfies, in spherical coordinates,

$\nabla\cdot\textbf{{P}}=-\rho\kappa_{\mathrm{R}}\frac{\textbf{{F}}}{c},$ (31)

where F is the vector flux. Solving this equation for the radiation energy
density $E_{r}$ requires a closure relation. We use the analytical $M_{1}$
closure in this model, since it is both purely local (removing the need for
any iteration) and exact for a single source. This closure expresses $P_{rr}$,
the radial component of P, in terms of $E_{r}$ as $P_{rr}=f_{rr}E_{r}$, where
$f_{rr}$ is the radial-diagonal component of the Eddington tensor. Under the
$M_{1}$ approximation, this is related to the reduced flux $\zeta=F/(cE_{r})$
as

$\zeta=\frac{5-2\sqrt{4-3f_{rr}^{2}}}{3}.$ (32)

The non-radial components of P that arise when applying the divergence
operator can be written in terms of $P_{\mathrm{rr}}$ and $E_{r}$ by using the
fact that the trace of P is equal to $E_{r}$. Subsequent substitutions allow
us to write an ODE for $\zeta$ given by

$\displaystyle\partial_{r}\zeta=$
$\displaystyle\frac{3\zeta\sqrt{4-3\zeta^{2}}}{5\sqrt{4-3\zeta^{2}}-8}\left[\frac{\partial_{r}\ln
F}{3}\left(5-2\sqrt{4-3\zeta^{2}}\right)\right.$ (33)
$\displaystyle\left.+\frac{2}{r}\left(2-\sqrt{4-3\zeta^{2}}\right)+\rho\kappa_{\mathrm{R}}\zeta\right].$

This ODE passes through a critical point when
$\zeta_{\mathrm{crit}}=2\sqrt{3}/5$, at a radial location $r_{\mathrm{crit}}$
that can be obtained numerically by solving

$3\partial_{r}\ln F+\frac{4}{r}+2\sqrt{3}\rho\kappa_{\mathrm{R}}=0.$ (34)

Our numerical procedure is to solve for the location of the critical point,
and then integrate outwards and inwards in radius.

Thus far our approach parallels that of Skinner & Ostriker (2013). Where we
diverge is that they treat the opacity $\kappa_{\mathrm{R}}$ as constant,
whereas for us the opacity is dependent on the radiation temperature
$T_{r}=(E_{r}/a_{\mathrm{R}})^{0.25}$, which itself can be rewritten in terms
of $\zeta$ (and $F$). This renders Equation (33) an ODE in one variable
($\zeta$), which can be integrated numerically to obtain $\zeta(r)$, and
subsequently $T_{\mathrm{r}}(r)$, and $\kappa_{\mathrm{R}}(r)$. We show the
resulting profiles for some examples of $\Sigma_{\mathrm{cloud}}$ and $\alpha$
in Figure 23.

Once we have obtained the profiles, it is straightforward to compute the
radiative force as $f_{\mathrm{rad}}=\rho\kappa_{\mathrm{R}}F/c$. The
corresponding forces of gravity on the gas are given by the sum of their
interaction with the source ($g_{\mathrm{*}}$) and self-gravitational
interactions ($g_{\mathrm{self}}$) as
$f_{\mathrm{grav}}=\rho(g_{\mathrm{sink}}+g_{\mathrm{self}})$. The former is
simply $g_{*}=GM_{*}/r^{2}$, and
$g_{\mathrm{self}}=G\int_{0}^{r}dM(r^{\prime})/r^{2}$, where
$\int_{0}^{r}dM(r^{\prime})=\int_{0}^{r}\rho(r^{\prime})4\pi{r^{\prime}}^{2}dr^{\prime}$
is the cumulative gas mass out to radius $r$. These can then be combined to
construct the Eddington ratio
$f_{\mathrm{Edd}}=f_{\mathrm{rad}}/f_{\mathrm{grav}}$. The numerical
implementation for this model is available at the following github repository:
https://github.com/shm-1996/ModelCloud.