Optical matter (OM) consists of arrays of micro or nanoparticles that are somehow bound and controlled by light (Burns et al., 1990). An OM able to self-assemble at will to develop solid technology is a long-standing goal in photonics (Parker et al., 2020). The background for OM is the particle manipulation by optical forces; the first results were applied to microparticles due to the lack of technology and the presence of thermal noise for smaller systems (Ashkin and Dziedzic, 1975; Ashkin and Dziedzic, 1989). Generally, OM comprises multiple particles subjected to electromagnetic forces that come from their mutually scattered light. However, this multiple scattering phenomena can be very complex and bring about unusual effects such as “non-reciprocal” forces, torque opposite to the illumination angular momentum, and non-conservative forces (Albaladejo et al., 2009; Parker et al., 2020). Appropriate control of OM by forces and torques could lead to programmable materials for optomechanical, rheological, and biological applications. In this respect, many works studied the optical binding between nanoparticles as a primary tool to develop OM. Some approached specific combinations of optical beams like those with programmable phase (Han et al., 2018; Nan et al., 2022). For example, light-induced rotation of objects holds potential for various applications such as sensing, cargo transportation, drug delivery, and micro/nanosurgery (Arita et al., 2020; Blázquez-Castro et al., 2020; Xin et al., 2020). Optical traps use the combination of beams as a potent characterization tool for material science and biophysics, as in Ref. (Erdogan et al., 2022), which uses electrostatic focusing to obtain the mass spectrum of SARS-CoV-2 and BoHV-1 virions. However, the high intensity at the focal spot may introduce laser heating, which is an issue for bio applications (Nan et al., 2022).

On the other hand, recent advances in THz technology call for new devices and materials that exhibit a non-reciprocal behavior for photonic networks and optical information processing (Shui et al., 2022). Non-reciprocal devices are a crucial component of modern communication technology. They are nowadays required for miniaturized electronic and photonic devices (Nature, 2020). One way to create optical non-reciprocity at the THz range is using magneto-optical (MO) systems like graphene, hexaferrites, and semiconductors (Kimel et al., 2022). For instance, Ref (Chochol et al., 2016). presented MO measurement of several samples of InSb with different carriers and carrier concentrations for low external magnetic field and room temperatures. With these advantages, dimers and trimers of InSb particles have been studied to enhance THz spectroscopy by forming electric and magnetic hotspots in the gap between them (Bakker et al., 2015; Sadrara and Miri, 2019). Anisotropic materials like InSb in OM would make it strongly dependent on the beam combinations, allowing for countless possibilities (Roichman et al., 2008). Recently, several efforts have focused on MO nanoparticle systems to shape OM and optical traps with a reasonable degree of control and accuracy, besides other relevant applications (Emile et al., 1992; Hu et al., 2013; Abraham Ekeroth et al., 2017; Edelstein et al., 2019; Edelstein et al., 2022). On a broader sense, magnetoplasmonics relates the plasmonic behavior of nanoparticles with the presence of external magnetic fields. The modulation and tunability of plasmonic resonances offered by magnetoplasmonics results auspicious for ultra-sensitive sensors and active plasmonic devices (Lacroix et al., 2012). In particular, the formation of stable optical binding between two small magnetoplasmonic particles has been lately studied (Edelstein et al., 2021). Equilibrium binding distances were predicted and found tunable by the incoming wave’s polarization state and the magnetic field’s magnitude. However, the model developed in that report is valid only for relatively small magnetic field values. Moreover, it predicts stable dimers only using alternating magnetic static fields and polarization angles to remove azimuthal, unbalanced forces. More importantly, that work needs to discuss possible rotations of the particles due to angular momentum transfer in the multiple scattering scheme.

In this paper, we study the formation of stable MO dimers for small nanoparticles in a complete framework involving all the possible optomechanical inductions. The dimer’s isotropic and anisotropic responses are assessed as a base for OM designs, i.e., under the presence/absence of an external magnetic field. This field can be of arbitrary magnitude in our model. The illumination consists of two counter-propagating plane waves with circular polarization, which simulates a simple optical trap in the vacuum. We found several possibilities to create stable dimers when the magnetic field is on. The beams do not exert net forces for reciprocal dimers but may exert torques on them. Nonetheless, there is a net radiation pressure and spin for the whole system when the static field is on, allowing complete control of the system’s movement. The results will enable one to infer that the mechanical variables can be used as near-field observables to explore the content of unknown samples. Conversely, they can be used to accurately control the dimer’s creation/destruction and its mobility. Finally, the spins predicted are explained in terms of the energy flows around the dimer, which constitutes a novel scattering-force effect for interacting particle arrays.

In the following, we assume two equal particles of the same non-reciprocal material immersed in the vacuum. In particular, we simulate nanoparticles made of n-doped Indium antimonide (n-InSb) (Palik et al., 1976) when a static magnetic field is along the z-direction, B = B z ̌ , see Figure 1. In this case, the permittivity tensor of InSb adopts the following form ϵ ̂ B = ϵ 1 B − i ϵ 2 B 0 i ϵ 2 B ϵ 1 B 0 0 0 ϵ 3 , (1) where ϵ 1 B = ϵ ∞ 1 + ω L 2 − ω T 2 ω T 2 − ω 2 − i Γ ω + ω p 2 ω + i γ ω ω c 2 − ω + i γ 2 , ϵ 2 B = ϵ ∞ ω p 2 ω c ω ω + i γ 2 − ω c 2 , ϵ 3 = ϵ ∞ 1 + ω L 2 − ω T 2 ω T 2 − ω 2 − i Γ ω − ω p 2 ω ω + i γ . (2) Here, ϵ _∞ is the high-frequency dielectric constant, ω _L is the longitudinal optical phonon frequency, ω _T is the transverse optical phonon frequency, ω p 2 = n e 2 / ( m * ϵ 0 ϵ ∞ ) defines the plasma frequency of free carriers of density n and effective mass m*, Γ is the phonon damping constant, and γ is the free-carrier damping constant. Noteworthy, n relates to the doping concentration in n-InSb; any modification of this doping would modify the spectral resonances of the system through the influence of ω _p on the dielectric tensor given by Eq. 1. The magnetic field enters in these expressions via the cyclotron frequency ω _c = eB/m*. In all calculations below, ϵ _∞ = 15.7, ω _L = 3.62 × 10¹³ rad/s, ω _T = 3.39 × 10¹³ rad/s, Γ = 5.65 × 10¹¹ rad/s, γ = 3.39 × 10¹² rad/s, n = 1.07 × 10¹⁷ cm⁻³, m*/m = 0.022, and ω _p = 3.14 × 10¹³ rad/s. Indium antimonide (InSb) is one example of the most widely studied polar semiconductors for magnetoplasmonic applications because it can be easily doped for sizable magnetic-induced effects (Moncada-Villa et al., 2015; Chochol et al., 2016). As reviewed in Refs. (Palik et al., 1976; Moncada-Villa et al., 2015), n-InSb is an exciting material that has two kinds of surface resonances in the absence of static field, namely, the phonon polariton (SPhP, higher-energy) and the plasmon polariton (SPP, lower-energy). For the reader’s suitability, Figure 2 shows the principal non-zero components of the permittivity tensor 1 as a function of the wavelength when B = 1 T. The surface resonances can be tracked from these spectra at the wavelengths that fulfill the corresponding Frölich condition (Abraham-Ekeroth, 2022).

For two small particles compared to the wavelength, we can make use of the approximation of discrete dipoles (Discrete Dipole Approximation, DDA) (Draine and Flatau, 1994; de Sousa et al., 2016). The usual DDA approximates a continuum system by a finite array of polarizable cells, thus constituting a non-relativistic integral solution of Maxwell equations. The cells acquire dipole moments in response to the local fields. These dipoles interact with one another via the appropriate Green functions (Novotny and Hecht, 2006). Then, proper cell polarizabilities can consider multiple-scattering scenarios. The principal limitation of the DDA concerns the handling of particle boundaries and interparticle distances since the geometry’s representation of the array has a minimum length scale equal to the inter-dipole spacing (Draine and Flatau, 1994; Yurkin and Hoekstra, 2007). Following the general criteria (Draine and Flatau, 1994), we represent the particles by dipoles, provided they are much smaller than the wavelengths under study. Similarly, we deal with inter-particle gaps relatively large enough to model a dimer with two dipoles. In such a case, the DDA simplifies considerably to p 1 = ϵ 0 α ̂ E 0,1 + k 0 2 α ̂ G ̂ 12 p 2 , (3) p 2 = ϵ 0 α ̂ E 0,2 + k 0 2 α ̂ G ̂ 21 p 1 , (4) where α ̂ is the polarizability tensor representing the particles. Thus this model assumes a linear relation between the dipole moment and field beside the approximation of small particles. The following definition automatically includes the radiative corrections necessary to fulfil the optical theorem (Ekeroth, 2019) α ̂ = α ̂ 0 − 1 − i k 0 3 6 π I ̂ − 1 (5) where α ̂ 0 is the so-called quasistatic polarizability, which can be given by α ̂ 0 − 1 = 1 V L ̂ + ϵ ̂ r − I ̂ − 1 (6) being V the particles’ volume, L ̂ = I ̂ / 3 is the electrostatic depolarization tensor specified for spheres or cubes, and ϵ ̂ r is the relative dielectric tensor. In the current work, we consider spherical nanoparticles. Note that the model omits magnetic polarizabilities since n-InSb results non-magnetic in the optical range of interest. G ̂ is the Green tensor defined as G ̂ r , r ′ = e i k 0 R 4 π R 1 + i k 0 R − 1 k 0 2 R 2 1 ̂ + 3 − 3 i k 0 R − k 0 2 R 2 k 0 2 R 2 R ⊗ R R 2 , (7) where R = r − r′, R = |r − r′|, and ⊗ denotes the exterior product (Novotny and Hecht, 2006). In particular, G ̂ 12 = G ̂ 21 in Eq. 7 refers to G ̂ evaluated at the particles’ centers. The system of equations 3 can be solved straightforwardly, leading to p 1 = ϵ 0 F ̂ E 0,1 + k 0 2 α ̂ G ̂ 12 α ̂ E 0,2 , (8) p 2 = ϵ 0 F ̂ E 0,2 + k 0 2 α ̂ G ̂ 12 α ̂ E 0,1 , (9) where we define F ̂ = α ̂ − k 0 4 G ̂ α ̂ G ̂ − 1 . In this work, a counter-propagating configuration is assumed as a superposition of two left-handed circularly polarized (LCP) plane waves with the same intensity I ₀ (Cameron et al., 2014; Edelstein et al., 2019), see Figure 1, namely, E 0 = E 0 2 x ̌ + i y ̌ e i k 0 z + x ̌ − i y ̌ e − i k 0 z . (10) This field is used in Eq. 8 to calculate the incident field at the particles’ positions E _0,1 and E _0,2.

The absorption cross section of the system can be calculated once the dipole moments are known by σ abs = k 0 ϵ 0 w E tot I m p 1 ⋅ α ̂ 0 − 1 p 1 * + p 2 ⋅ α ̂ 0 − 1 p 2 * (11) where w E tot = ϵ 0 | E 0 | 2 . The i-component of the forces exerted on each particle can be obtained from the time-averaged force within the Rayleigh approximation (Chaumet and Nieto-Vesperinas, 2000). This is F 1 , i = 1 2 R e p 1 t ∂ i E * r , ω | r = r 1 (12) F 2 , i = 1 2 R e p 2 t ∂ i E * r , ω | r = r 2 (13) where the derivatives of the total field ∂ i E ( r , ω ) | r = r n at the dipoles’ positions r _n , n = {1, 2} can be obtained from (Chaumet and Billaudeau, 2007): ∂ i E r , ω | r = r 1 = ∂ i E 0 r , ω | r = r 1 + + k 0 2 ϵ 0 ∂ i G r , r 2 r = r 1 p 2 (14) ∂ i E r , ω | r = r 2 = ∂ i E 0 r , ω | r = r 2 + + k 0 2 ϵ 0 ∂ i G r 1 , r r = r 2 p 1 (15) The total force exerted on the dimer results from adding the force components for each particle, namely, F _tot,i = F _1,i + F _2,i . In particular, the net radiation pressure for the dimer under the illumination given by Eq. 10 is defined by taking i = 3, or the z components, as F tot,z = F 1 , z + F 2 , z (16) Another useful mechanical variable is the binding force, which in the present case is defined as Δ = F 1 − F 2 ⋅ n ̌ (17) where n ̌ = r 2 − r 1 | r 2 − r 1 | is the dimer’s versor. For perpendicular configuration, Figure 1B, it is convenient to define the following tangential force for the system τ = F 1 ⋅ ϕ ̌ 1 + F 2 ⋅ ϕ ̌ 2 (18) where ϕ ̌ n correspond to the tangential versors of the variable ϕ for each particle. The optical torques can also be calculated, as given in Ref. (Chaumet and Rahmani, 2009): N spin,1 = 1 2 ϵ 0 R e p 1 × α ̂ 0 − 1 * p 1 * (19) N orb,1 = r 1 × F 1 (20) N 1 = N spin,1 + N orb,1 (21) N spin,2 = 1 2 ϵ 0 R e p 2 × α ̂ 0 − 1 * p 2 * (22) N orb,2 = r 2 × F 2 (23) N 2 = N spin,2 + N orb,2 (24) The definitions of the orbital and spin torques were discussed previously in Refs. (Chaumet and Rahmani, 2009; Nieto-Vesperinas, 2015a; Nieto-Vesperinas, 2015b), among others. The spin torques are always defined with respect to the centers of the particles. Otherwise, the reference system is located at the dimer’s center of mass, and orbital torques are set. Thus, the total torque exerted on the dimer is N orb = N orb,1 + N orb,2 (25) N spin = N spin,1 + N spin,2 (26) N tot = N 1 + N 2 = N orb + N spin (27)

Since we are interested in the near-field interactions between the particles, the study focuses on an example for which the interparticle’s gap equals one particle diameter (g = 2R _p ); see Figure 1. We add complementary examples for other values of g in the Supplementary Material (SM).

In this section, all the optical variables were scaled by the proper factors to make them dimensionless. The following characteristic magnitudes, namely, w E tot , A p = π R p 2 , V p = 4 3 π R p 3 , and V int = 4 3 π | r 2 − r 1 | 3 , where R _p is the particles’ radius, redefine the variables as Q abs = σ abs 2 A p for the absorption efficiency, F rad = F tot,z w E tot A p and Δ ′ = Δ w E tot A p for the radiation pressure and the binding forces, τ ′ = τ w E tot A p for the tangential force, and N spin ′ = N spin w E tot V p for the spin torque. An example only shows the variable N orb ′ = N orb w E tot V int for the orbital torque in the SM since it gave negligible results compared to the spin N spin ′ unless g is minimal, see Supplementary Figures S1, S2 for details. We calculate the scaled Poynting vector as S = 1 2 I 0 R e { E × H * } where the magnetic field H comes from an DDA equation similar to that for E (Novotny and Hecht, 2006). The curl of S is calculated using an appropriate three-dimensional mesh around the system’s near-field. We define a “stable” dimer whenever the set of conditions Δ ′ = 0 κ x ′ = d Δ ′ d x < 0 τ ′ = 0 κ ϕ ′ = d τ ′ d ϕ < 0 F rad = 0 (28) are approximately fulfilled in our model (Dholakia and Zemánek, 2010). Here, x = k ₀ R and we define two more dimensionless numbers, i.e., κ x ′ and κ ϕ ′ , which mean the radial and azimuthal “spring” constants. The last condition is not strictly necessary for dimer’s creation and can be relaxed. Although the movements of the particles are not approached because their velocities and accelerations are considered zero, this model constitutes a useful approximation to approach the creation and dimers’ dynamics (Lamothe et al., 2007; Haefner et al., 2009; Miljković et al., 2010; Liaw et al., 2016).

By a simple dipolar model, this work explores the behavior of small nanoparticle dimers when magneto-optical materials like n-doped InSb and moderate magnetic fields are used. Two counter-propagating waves with equal circular polarization are used as illumination to simulate a simple optical trap with neither net gradient nor scattering forces. Our results show that the system can be thoroughly characterized by observing its mechanical inductions, provided these latter depend on the near field. Besides, we found no stable dimers when the dimer is aligned with the illumination since the inter-particle force only leads to repulsion. On the contrary, under “perpendicular” alignment and in the presence of a magnetic field, we found dimer’s stability is possible, at least under this “static” model for which the particles’ velocities and accelerations are not considered. As the results strongly depend on the magnetic field’s presence, this study constitutes a novel background to build OM or photonic molecule nano factories and control their movements, or conversely, to study samples containing this class of multimers.

Finally, we give an original explanation for the appearance of particle spins based on the energy flows. This interpretation offers satisfactory results from the “scattering” forces produced by the interaction between the particles. Gradient forces were also investigated but showed no appreciable influence on the spins’ appearances (results not shown in this work).

It is worth mentioning that there is a compromise between the particle size and inter-particle gap to maintain the maps in Figures 4, 5 “invariant” or scalable. In the region between R _p = 0.05 and 0.65 μm for a condition g = 2R _p , or when R _p = 0.2 μm for gaps above g = 0.1R _p , very similar maps to those reported in Figures 4, 5 are obtained. There are only variations in the relative colormap intensities. As particle sizes increase compared to g, some hotspots that appear at special values of ϕ start to accentuate or deform. The shapes start distorting more from particle radii R _p ≅0.7 μm and beyond, or for interparticle gaps shorter than g = 0.1R _p when R _p ∼ 0.2 μm. Remarkably, the maps for Δ′ may provide newer “zero” lines according to the distortions. However, when one increases g and fixes the particles’ size, the original shapes are recovered. Thus, we can state a “rule of thumb”: our results apply to nanoparticle dimers with particle sizes between 0.1 and 1.4 μm as long as g remains in the order or larger than R _p . For other lengths, the model presented in Section 2 still works but the maps would change their shape. Nonetheless, smaller particles or tiny gaps may require a size or non-local correction to the dielectric tensor of the system to recover realistic results (García de Abajo, 2008; Scaffardi et al., 2013). This work provides a few examples to illustrate the general model. Exhaustive studies about the parametric variations of geometry or particle constitution are considered beyond the scope of this paper.