Abstract

Controlling individual atoms using lasers 1 , ion traps 2 and scanning probe tips 3 has transformed our understanding of matter and enabled breakthroughs in quantum science 4 , 5 , 6 . Extending this control into three-dimensional (3D) solids and across mesoscopic scales, however, remains a foundational challenge. Electron irradiation in electron microscopes is known to induce atomic displacements 7 , and atomic manipulation has been proposed 8 and demonstrated 9 , 10 . Yet repeated and deterministic control has remained elusive 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 . Here we demonstrate deterministic atomic engineering in a 3D crystal, creating ordered arrangements of more than 40,000 user-defined defects within minutes across a 150 nm × 100 nm × 13 nm volume. By steering individual Cr atoms in the magnetic semiconductor CrSBr into selected interstitial sites using an electron beam directed with sub-20-pm-scale accuracy, we create vacancy–interstitial complexes. The resulting impurity array forms a mesoscale crystal embedded within the host lattice, a new form of engineered artificial matter that remains stable at room temperature and outside the microscope. By tracking Cr atom displacements, we identify conditions under which the defect structures are predictable. Our calculations suggest that these defects form correlated impurity states with intra-defect optical transitions and inter-defect kinetic and Coulomb interactions. This establishes a generalizable platform for atomic defect engineering at mesoscopic, and potentially macroscopic, scales, opening opportunities for scalable quantum technologies, including deterministic colour-centre placement, quantum simulation of many-body lattice models and atomic-scale manufacturing.

Main

Atomic-scale manipulation has emerged as a defining ability of modern quantum science and materials research and has profoundly shaped the fields of quantum simulation 18 , 19 , computing 20 and high-density storage 21 . Central to these achievements are techniques that use optical 1 and oscillating electric fields 2 to trap atoms in vacuum, or that use scanning probe tips to reposition atoms on surfaces 3 . Extending this level of control into three-dimensional (3D) solids and across mesoscopic scales remains a foundational challenge that has deep implications for designing artificial quantum materials and for scaling quantum technologies.

Ions or electrons offer strategies for structural modification to create functional defects in 3D solids. Although colour centres 22 or other atomic-scale spin–photon interfaces can be formed this way, scalable integration into photonic architectures 23 requires a higher degree of control of the position and nature of the defects produced than is available at present. Direct atomic engineering in 3D solids is theoretically conceivable using the sub-ångström focused electron beam in a scanning transmission electron microscope (STEM). Manual displacement of dopants and atomic rearrangements have been achieved in both thin crystals (such as graphene) 8 , 9 , 10 , 11 , 13 , 14 , 16 , 17 and thicker materials 12 , 15 , 24 , 25 . A broader vision is outlined in a commentary article 26 , with previous work and challenges summarized in recent reviews 27 , 28 . However, scalable and repeatable atomic engineering with deterministic control and spatial programmability has yet to be demonstrated. A technique that can engineer optically active spin and charge defects and artificial lattices with controlled position, and hence programmable spacing and tunable interactions, would allow us to imagine artificial structures with tailored symmetry that transcend the constraints of naturally occurring atomic arrangements. This ability could establish solid-state analogues to cold-atom quantum simulators that allow access to multi-orbital many-body Hamiltonians with strong Coulomb and exchange interactions, regimes currently beyond the reach of both experiment and conventional theory.

Here, we demonstrate deterministic atomic engineering of distinct, spatially programmable defect complexes throughout a 3D volume of a crystal. We perform these structural modifications by directing the movement of Cr atoms within CrSBr, an air-stable magnetic semiconductor. We use a focused electron beam to engineer user-defined defect types at precise locations across mesoscopic length scales, creating tens of thousands of defects by directing the beam rapidly and with few-picometre-scale precision. We show how vacancy–interstitial defects of different types are selectively generated, and we identify the underlying displacement mechanism and quantify and discuss the necessity for few-picometre-precision electron delivery. We perform first-principles calculations to predict that these defect complexes should possess electronic structures analogous to the NV − centre in diamond. We finally create defect arrays with different spacings and calculate the collective behaviour that the resulting ensembles are expected to host. The results indicate that these arrays may offer a promising platform for exploring many-body physics in the solid state.

Deterministic 3D atomic engineering

Figure 1 a shows a high-angle annular dark field (HAADF)-STEM image of 24 precisely written columns of engineered defects in a 16-layer (12.7 nm-thick) CrSBr crystal. Each column contains one vacancy–interstitial pair per layer, constituting a total of >300 defect complexes arranged in a 3D volume. As we discuss below, these complexes are electronically decoupled, with finite Coulomb interaction but negligible out-of-plane hopping.

a, HAADF-STEM image (unprocessed data) of a 4 × 6 grid of V Cr +Cr int delta 2 electron-beam-created defect complexes. In the unaltered lattice, the S/Br columns are bright, and the Cr columns have lower intensity. The box outlines one beam-generated defect column. The total time per defect column is about 1.25 s. b, HAADF-STEM image of a 6 × 2 grid of V Cr +Cr int delta 1 defect complexes. c, Corresponding schematics and images showing the effect of a periodic 1D beam motion with a frequency of 0.5 kHz for 1 s in the ±b- and ±a-directions (left and right, respectively). Cr atoms are displaced into the van der Waals gap creating a Cr interstitial Cr int and a leftover Cr vacancy V Cr at each layer of the crystal. d, Autonomous atomic engineering at scale by moving the electron beam between target locations. At each target location, the beam is positioned with few-picometre precision; then a beam action is performed to steer a column of Cr atoms into the chosen interstitial locations. This forms a column of vacancy–interstitial pairs, with each pair acting as an impurity state. The nature of the defect column can be different at each location (bottom). e, HAADF-STEM image of an engineered star. All patterns are created and imaged with a beam energy of 200 keV and a beam current of 20 pA. Scale bar, 2 nm (a, b, e).

CrSBr offers an excellent model system for these experiments. It is an air-stable 29 layered magnetic semiconductor 30 with weak interlayer interactions 31 , as seen through its quasi-1D electronic structure, its monolayer-like optical properties 32 and its sharp defect emission even in bulk 33 . Most importantly, CrSBr is known to exhibit electron-beam-induced Cr mobility. At 200 keV, Cr atoms can move between adjacent lattice and interstitial sites in both multilayer 34 and bilayer 35 CrSBr under irradiation. The nature of the defects and structural rearrangements in CrSBr can be identified readily because in plan view, the lattice exhibits two distinct atomic-column types containing either only Cr atoms or both S and Br atoms.

The patterns in Fig. 1a were each created by a rapid, spatially programmed 1D beam trajectory that spans three atomic columns (S/Br–Cr–S/Br) in a back and forth motion along the b-axis (Fig. 1c). The contrast in the image indicates that this motion has displaced Cr atoms from their initial locations (δ0) into adjacent S/Br column locations (δ2), forming a vertical stack of complexes we label V Cr +Cr int delta 2, with the interstitial presumably within the van der Waals gap, as discussed below. The Cr displacement direction follows the beam path: rotating the 1D beam trajectory by 90 ° along the a-axis generates a different complex (Cr in location δ1) which we label V Cr +Cr int delta 1 (Fig. 1b).

Achieving this control required a specific, automated atomic engineering work flow that is described in Fig. 1d, Supplementary Information sections 1–2 and Supplementary Figs. 1–6. The key attributes are dose-efficient, few-picometre-precision beam positioning, rapid beam motions and pre-programmed movements between target sites. These abilities allow the rapid creation of complex and varied 3D defect structures with lattice-level precision (Fig. 1e). As we discuss below, reproducibility through the thickness of the crystal is achieved through electron channelling, enabling consistent formation of defects in 3D arrays (Supplementary Information section 3 and Supplementary Figs. 7–11).

Control of defect type

Figure 2 shows the structural outcomes of 10 distinct beam trajectories, each producing a unique defect configuration that is expected to display its own electronic character. For each beam trajectory type, 30 experiments were carried out, a HAADF-STEM image of each area was acquired immediately after beam engineering, and these images were averaged. The trajectories consisted of periodic 1D motions along three directions: vertical along the crystallographic a-direction (±a, 0), horizontal along the b-direction (0, ±b) and diagonal (±a, ±b). For vertical and horizontal paths, the beam is moved across either two (Cr–S/Br) or three (S/Br–Cr–S/Br) atomic columns. Diagonal motions are performed across two adjacent Cr columns (Cr–Cr). It is clear from the image contrast that the horizontal and vertical motions have steered Cr atoms into neighbouring S/Br columns along the beam direction, as indicated by increased HAADF intensity at these sites. Two-column motions enable selective placement into any of the four nearby S/Br columns, whereas three-column motions induce displacement into two adjacent sites. Diagonal motions result in Cr atoms occupying opposing S/Br columns orthogonal to the beam motion axis.

HAADF-STEM images collected after performing electron beam motions along the vertical, horizontal and diagonal directions. Each image represents an average of 30 individual experiments. For each experiment, the beam is moved for 1 s at a frequency of 0.5 kHz with a beam energy of 200 keV and a beam current of 20 pA. The CrSBr thickness is 16L (12.7 nm). Scale bar, 0.5 nm.

Importance of picometre positioning

Achieving this degree of control over Cr displacement required the sub-20-pm positioning accuracy made possible by our beam-control algorithms. We assess the importance of accurate electron beam positioning by performing 5,694 of the ±b, 1 s beam motions, each randomly positioned on the sample and followed by a small-field-of-view HAADF-STEM image at each location (Fig. 3a, Supplementary Information section 4 and Supplementary Figs. 12–14). By relating the image to the position of the beam motion within the Wigner–Seitz cell in the Cr sub-lattice, we can visualize the cause–effect relationship between beam position and the atomic structural response, with 10 pm resolution (Fig. 3b). Even offsets of <15 pm induce measurable intensity changes in neighbouring atomic columns (Fig. 3c,d). For example, diagonal offsets from the target centre Cr column (Cr {c} ) displace Cr atoms in the corresponding diagonal columns (Cr {d1−d4} ), whereas ±a-direction offsets shift Cr atoms towards the S/Br columns (S/Br {u , d} ). Gaussian fits to the spatial intensity distribution of Cr {c} indicate that <14.1 pm precision is required to maintain 98th-percentile intensity changes (Fig. 3d). We suggest that in reality positioning requirements are even more stringent, considering effects on adjacent atom columns. Although electron channelling helps to guide electrons, precise algorithmic beam control to attain this level of precision is crucial for defect reproducibility.

a, Wigner–Seitz cell in the Cr sub-lattice with 5,694 random electron beam positions, each representing the centre of a 1 s long ±b beam motion. b, Corresponding mosaic of averaged HAADF-STEM images, tiled at a spatial resolution of 10 pm × 10 pm. c, Integrated atomic-column intensities for the central Cr site and its nearest-neighbour diagonal Cr {d1−d4} and S/Br {l, r, t, d} columns. d, Two-dimensional Gaussian fit to the intensity distribution of the central Cr {c} atom column, along with horizontal and vertical cross-sections. A 2% intensity threshold corresponds to required positioning accuracies of σ x = 13.1 pm and σ y = 14.1 pm.

Defect-induced impurity states

The ability to select the defect type at each location as shown in Fig. 2, coupled with the full programmability of defect superlattice spacing and pattern, with examples shown in Fig. 4a, suggests that it is possible to engineer structures having sites with different local electronic properties and systematically explore collective electronic phenomena that could emerge from patterned defect architectures that are non-locally coupled. As one of many intriguing and representative examples, we focus on the V Cr +Cr int delta 2 defect complex. To investigate the properties of the individual defects, we show in Fig. 4b the electronic structure of the V Cr +Cr int delta 2 defect complex calculated within quasiparticle self-consistent GW (QSGhatW) theory 36 , 37 , 38 . More details for this and other defects are provided in Supplementary Information section 5 and Supplementary Figs. 15–17. The defect introduces five flat impurity states (Fig. 4b, black bands and arrows), three fully occupied and two empty. Three of them originate from three atomic-like orbitals: two p x -like (dangling-bond-like and rotated) states on neighbouring S atoms adjacent to the Cr vacancy (p x (1), p x (2)), and a d x y -like orbital on the interstitial Cr atom (Fig. 4c), reminiscent of the NV − centre. Here, the two S p x (1,2) states form a bonding–anti-bonding pair p x ^(b/a)=frac1sqrt2(p x (1) +/- p x (2) ). The anti-bonding orbital further hybridizes with the Cr d x y state, such that the lowest of these three impurity states is of pure S p x ^(b) character, whereas the two upper impurity states are of mixed S p x ^(a) and Cr d x y character (Supplementary Information section 5 and Supplementary Fig. 18).

a, HAADF-STEM images of superlattices of V Cr +Cr int delta 2 defect complexes with varying periodicity. The insets show the corresponding FFTs. b, QSGhatW electronic band structure of V Cr +Cr int delta 2 with host bands (grey) and five defect bands (black). The k-space contribution of one of the five excitonic bound states (X_D 2) is highlighted (cyan). c, Orbital and electronic level structure of V Cr +Cr int ^delta 2. d, Distance-dependent intra- and inter-defect Coulomb and hopping energies. The vertical line indicates the minimal inter-defect separation (the unit cell distance). Grey shading is a guide to the eye. e, Bound exciton-level diagram from BSE calculations. Scale bar, 2 nm (a).

To assess the possibility of collective phenomena when these defects are placed in arrays, we estimate intra- and inter-defect electron hopping t and Coulomb interaction U parameters within the basis of the three impurity orbitals p x (1) ,p x (2) and d x y (Fig. 4d and Supplementary Table 1). We find that the intra-defect hopping matrix elements between the three impurity orbitals range from 80 meV to 180 meV, and the Coulomb interaction strengths range between 600 meV and 1,000 meV for intra-orbital terms and are around 250 meV for inter-orbital interactions. These values point to strong electronic correlations within this defect. The partial occupation of the mid-gap impurity states suggests the formation of tightly bound intra-defect excitons, which is confirmed by full ab initio solutions to the Bethe–Salpeter equation (BSE) based on QSGhatW inputs. The resulting excitonic resonances appear in the near- to mid-infrared range (Fig. 4d, Supplementary Information section 5 and Supplementary Figs. 15–17), among which the X_D 2 feature corresponds to a strongly localized exciton, with its wavefunction projected onto the band structure in Fig. 4b (cyan). The isolated V Cr +Cr int delta 2 defect thus constitutes a correlated quantum object.

For the inter-defect Coulomb matrix elements, we estimate a r −1 -like dependence (Fig. 4d) with interaction strengths of several 100 meV. Inter-defect hopping parameters t remain above 10 meV for separations <1.4 nm, indicating coupling between neighbouring defects. Similarly to magneto-excitons in CrSBr 39 , interlayer impurity hopping is suppressed in the antiferromagnetic (AFM) phase, effectively confining hopping to 2D. In the FM phase, spin alignment enables 3D hopping. This is an exciting finding, as it allows the defect lattice physics to be described within a strongly correlated low-energy Hamiltonian of the form

with i and j orbital indices on the same defect as well as between different defects. t i j are the intra- and inter-defect hopping matrix elements, U i i are the defect and orbital local Coulomb interactions and U i j are the inter-orbital and/or inter-defect non-local Coulomb interactions. Depending on the orbital character and inter-defect spacing, U i i /t exceeds 10, placing the system deep in the strong-coupling regime (U i i ≫t). This interaction strength is relevant to engineering correlated quantum phases such as Mott insulating behaviour 40 , spin–charge separation and exotic orbital ordering in anisotropic multi-orbital lattices 41 , 42 , 43 .

Three-dimensional mesoscopic defect crystal

In building systems to probe correlated physics, it is essential that the atomic engineering technique has high scalability. In Fig. 5a,b, we demonstrate this property for our CrSBr patterning with an array of V Cr +Cr int delta 2 complexes at 4,582 (79 × 58) programmed positions, resulting in more than 40,000 individual defects within a 150 nm × 100 nm × 13 nm volume, written in about 40 min. Additional patterns with tighter spacings and other user-selected defects are shown in Supplementary Information section 6 and Supplementary Fig. 19. Sharp superlattice Bragg reflections in nanobeam electron diffraction (Fig. 5c) show that the written defects form a periodically modulated crystal-scale superstructure (Supplementary Information section 7 and Supplementary Fig. 20). This is a hallmark of mesoscale order and underpins the generation of artificial matter.

a, b, Grid of 79 × 58 locations in a volume of 150 nm × 100 nm × 13 nm. The defect spacing is four Cr–Cr distances in the b-direction and five in the a-direction. The total write time was about 40 min. The CrSBr thickness is 16L (12.7 nm). c, Diffraction pattern with a 50-nm diameter electron beam shows sharp superlattice Bragg reflections (inset) revealing the periodically modulated crystal-scale superstructure, a hallmark of mesoscale order. The first-order superlattice reciprocal vectors are given by bfg 1 =(bfa * /10) and bfg 2 =(bfb * /8) reflecting the total number of atom columns between defects. The pattern is created with a beam energy of 200 keV and a beam current of 60 pA with 0.5 s per site. Scale bars, 0.1 μm (a); 2 nm (b).

Deterministically arranging thousands of strongly interacting, multi-orbital defects with expected Coulomb repulsion U i i /t > 10 suggests a programmable platform for simulating quantum many-body Hamiltonians at a scale beyond the reach of both current classical techniques and early-stage quantum computers 44 . Moreover, we note that these lattices remain structurally stable outside the microscope (Supplementary Information section 8 and Supplementary Figs. 21–23), enabled by the material thickness and the sub-surface nature of the atomic engineering.

Three-dimensional defect-formation dynamics

The remarkable efficiency of directed Cr displacement probably originates from the interplay between the knock-on threshold and octahedral bonding environment of Cr and the quasi-1D electronic structure of CrSBr. Although Cr atoms are strongly bonded 33 , CrSBr exhibits highly anisotropic conduction bands dominated by Cr d-orbitals, strong exciton–phonon coupling and a Peierls-like lattice instability 32 . Moreover, the ground-state exciton also has a dominant Cr d-orbital character, leading to excitonic charge density localized around Cr atoms 32 . Further, long-lived dark excitons exist in CrSBr 45 that can induce a non-equilibrium energy landscape for typical beam currents. Together, these factors soften the lattice, which, we suggest, reduces the effective energy barrier for Cr migration under irradiation.

To probe the defect-formation mechanism, we systematically study the Cr displacement dynamics as the vertical stack of defects forms at the irradiated location. Figure 6a shows the change in the structure as a function of irradiation time, in images recorded after irradiating with the oscillating ±b trajectory for times from 5 ms to 1 s. Cr column intensities decrease with exposure, whereas S/Br column intensities increase, indicating progressive displacement. Our STEM image simulations suggest, based on the measured column intensities, that for a 1 s beam motion time, 75 ± 25% of Cr atoms are no longer at their initial sites (Supplementary Information section 8 and Supplementary Fig. 24). The ability of the beam to displace atoms deep down in a 3D crystal can be explained by the channelling of the beam through the lattice, as described in Supplementary Information section 3 and Supplementary Figs. 7–11.

a, Averaged HAADF-STEM images from 30 experiments collected after a ±b beam motion at different total durations with a frequency of 0.5 kHz, beam energy 200 keV and beam current 20 pA. The CrSBr thickness is 16L (12.7 nm). b, Real-time monitoring of Cr atom movement using fast (10 μs) HAADF detector readout over 1.5 s. The detector output signal is shown in different colours when the beam is on the (S/Br) (1) (light red), (S/Br) (2) (dark red), and Cr (grey) atomic columns. Initial and final local HAADF images are shown. c, Atom streaking showing the time evolution of Cr atom movement between three atomic columns. d, DFT calculated migration energy path of Cr atom from δ0 to δ2. e, Top, average of 25 experiments at 200 keV showing time-dependent Cr atom displacement dynamics. Bottom, integrated intensity of Cr, (S/Br) (1) and (S/Br) (2) atomic columns. The inset shows three-level rate equation model. f, Time-dependence at 100 keV. g, The Cr displacement model across multiple layers, nucleating at the entrance of the electron beam and cascading through the bulk of the material. a.u., arbitrary units. Scale bar, 0.5 nm (a).

The Cr displacement dynamics can be probed at higher (2 ms) time resolution if we monitor the HAADF detector output continuously during the oscillating beam motion (Fig. 6b). We integrate the detector signal over all time intervals during which the beam interacts with each of the three atomic columns, with the benefit that this provides a spatiotemporal map of the configuration of all three atomic columns throughout the scan interval, rather than only atomic density at a single point. We refer to this measurement modality as ‘atom streaking’. This allows us to track the movements of atoms out of the central Cr column (position δ0) and on to the two neighbouring S/Br columns (positions ±δ2) (Fig. 6c). The movement of individual Cr atoms appears between Cr and both S/Br columns, consistent with a density functional theory (DFT)-calculated migration energy profile, which shows a local energy minimum of 672 meV for the δ2 interstitial site (Fig. 6d, Supplementary Information section 8 and Supplementary Figs. 22 and 23). Although the Cr atom stabilizes temporarily in this minimum, the beam evidently provides sufficient energy to overcome the barrier and drive Cr migration between lattice sites. This reversible motion is consistent with the interstitial-type displacements observed here and would not be expected for vacancy-assisted mechanisms. We further calculate that such vacancy-assisted cases are energetically not favourable and would result in strong lattice distortions, which we do not observe (Supplementary Information section 8 and Supplementary Figs. 25–27). Our calculation suggests that formation of the Cr interstitial takes place through a bond-switching mechanism: one of the original four Cr–S bonds is preserved during displacement, whereas the other three Cr–S bonds are replaced by three Cr–Br bonds, including one to the adjacent layer (Fig. 6d). This minimal bond disruption facilitates efficient Cr migration with a high defect stability.

The Cr migration dynamics become even more evident when we study ensemble behaviour, for which we average 25 long-duration (4 s) experiments (Fig. 6e,f). At 200 keV, the averaged signal displays a characteristic exponential decrease in Cr intensity, accompanied by an exponential increase in S/Br signal, as expected for a three-level system. (Fig. 6e, inset, Supplementary Information section 8 and Supplementary Fig. 28). This confirms an atom-number-conserving window of time during which Cr atoms are steered between sites with directional control. This displacement mechanism strongly depends on beam energy. At 100 keV, the Cr intensity decreases linearly (Fig. 6e); the atom-number-conserving window disappears at lower energies and directional control is lost (Supplementary Information section 9 and Supplementary Figs. 29 and 30). At high beam energies, we suggest that knock-on momentum transfer induces bond softening and enables directional displacement by bond exchange along symmetry-preferred axes, but at lower energies we propose that both knock-on activation and electronic selectivity are suppressed, resulting in slower, nondirectional motion. At high beam energies, the Cr displacement rate is relatively high, allowing Cr atoms to move between defined sites before damage occurs and producing the clear three-level behaviour. We suggest that the atom-number-conserving displacement mechanism persists at all beam energies we studied, but at lower keV the Cr intensity decreases linearly (Fig. 6e) and becomes too slow to outpace general lattice damage, making directional motion unresolvable and giving the appearance of material loss (Supplementary Information section 9 and Supplementary Fig. 30). Thus, the observed energy dependence reflects a rate threshold rather than a change in mechanism. Rate equation modelling at 200 keV suggests cross-sections in the sub-barn to barn regime consistent with near-threshold displacement pathways, supporting a mixed process in which inelastic electron excitation transiently lowers the activation barrier while sub-threshold elastic momentum transfer provides the directional impulse (Supplementary Information section 8 and Supplementary Fig. 28).

We finally note that the DFT simulations (Supplementary Information section 9, Supplementary Fig. 31 and Supplementary Table 2) suggest that after a Cr atom is displaced out of one layer (for example, along +b), the adjacent Cr atom in the layer above or below will shift preferentially in the opposite direction (for example, −b) to minimize interstitial energy (Fig. 6e). This result indicates a cascading displacement mechanism that propagates through the thickness of the crystal.

Conclusion

We have demonstrated atomic engineering at mesoscopic scales in few-nanometre-thick CrSBr. Our approach combines deterministic, automated and 3D defect engineering within crystalline solids and, therefore, complements ion implantation (scalable but not atomically precise), scanning tunnelling microscope lithography (surface-limited) and optical tweezer or ion-trap arrays (subject to detrapping).

Looking ahead, we anticipate that the level of control shown here offers a pathway to quantum simulations with tunable many-body Hamiltonians, including multi-orbital configurations, strong Coulomb interactions and arbitrarily structured lattices that extend beyond the abilities of cold atoms, ion traps or moiré systems. Deterministic defect placement could also address challenges in scalable quantum computing, dense magnetic memory, and atomic-scale logic devices. As the samples, with defect arrays on a suspended membrane, are compatible with optical, electrical and magnetic measurements, we propose that our approach establishes a foundation for future correlative studies of the optical, spin and quantum properties of engineered defect patterns.

Our beam control strategies are, in principle, compatible with diverse materials and geometries, and potentially scalable to atomic structures comprising up to millions of precisely positioned defects, provided that additional algorithmic beam-control and probe-correction strategies are implemented. In particular, out-of-plane drift, tilt correction over large sample areas and compensation for dynamic aberrations such as astigmatism and spherical aberrations will be key to maintaining beam lock-on over extended durations. Furthermore, optimization of beam parameters such as probe current, beam energy and the emitted spatiotemporal dose profile offers additional room for improving writing speed. We believe that these developments could enable deterministic atomic engineering across micrometre-scale areas in less than 1 day.

Although CrSBr serves as a model, similar beam-atom displacement and defect-formation mechanisms are expected in a broader class of materials. We expect these effects to arise in compounds that offer energetically accessible local transformations, particularly those with sites that can accommodate interstitial configurations or that undergo beam-induced phase transitions. These criteria are not limited to layered systems and extend to non-layered compounds with flexible bonding environments or metastable coordination geometries. In general, materials known to exhibit structural or electronic phase transformations under external stimuli are promising candidates for extending beam-driven atomic engineering.

Beyond quantum material design, these results offer high-precision measurement modalities in STEM for quantitative studies of electron beam–solid interactions 46 , 47 , 48 , 49 , 50 , 51 , 52 , 53 , using the strategies shown in Fig. 6. Overall, we suggest that these results lay the foundation for programmable matter, in which atoms are placed with intent, and functionality is engineered from the atom up.

Methods

Crystal growth

The crystals were fabricated using chemical vapour transport as described in ref. 34 . CrSBr single crystals were synthesized by direct reaction of chromium (99.99%, −60 mesh, Chemsavers), bromine (99.9999%, Sigma-Aldrich) and sulfur (99.9999%, Stanford Materials). Bromine and sulfur were used with a 4 % excess over the stoichiometric ratio, corresponding to a total target mass of 15 g of CrSBr. The elements were loaded into a quartz ampoule (35 mm diameter, 220 mm), and the charge was cooled with liquid nitrogen to minimize bromine evaporation. The ampoule was sealed using an oxygen–hydrogen welding torch. The ampoule was sealed under high vacuum (<1 mPa) using a diffusion oil pump and a liquid nitrogen cold trap for evacuation.

The sealed ampoule was gradually heated from room temperature to 700 ° C over 5 days in a crucible furnace, while the opposite end of the ampoule was maintained below 200 ° C. This heating cycle was repeated until all liquid bromine had reacted.

Crystal growth was performed in a horizontal two-zone furnace. Initially, the growth zone was maintained at 900 ° C and the source zone at 700 ° C for 2 days. Following this, the thermal gradient was reversed: the source zone was ramped to 900 ° C and the growth zone cooled to 800 ° C, and these conditions were held for an additional 10 days. The resulting CrSBr crystals, typically up to 5 × 20 mm in size, were extracted from the ampoule inside an argon-filled glovebox.

Sample fabrication

Samples for STEM imaging were fabricated by mechanical exfoliation of CrSBr crystals onto SiO 2 /Si substrates. Atomic force microscopy and optical phase contrast were used for layer-precise thickness verification. Selected flakes were then transferred to STEM grids using a dry polymer-assisted transfer using cellulose acetate butyrate (CAB). For this transfer, the chip was then subjected to spin-coating with a thin layer of CAB, a polymer chosen for its low glass transition temperature of approximately 130 ° C, which allowed for baking without damaging the CrSBr flakes. A CAB solution in ethyl acetate (24 g/100 ml) was used, and the chip was subjected to spin-coating at 1,500 rpm for 60 s with an acceleration of 500 rpm s −1 , yielding a polymer thickness of approximately 15 μm. The coated chip was baked at 80 ° C for 6 min under a petri dish to ensure uniform baking and to prevent the surface from drying out.

A scalpel was used to cut a small region around the desired flake under an optical microscope. This polymer stamp was transferred mechanically with tweezers onto a holey STEM grid (Norcada NH050A2 or NH005D05). The STEM grid was coated with a 5 nm thin layer of carbon before the transfer to reduce charging effects during STEM imaging and defect engineering. After the polymer stamp transfer, the chip was baked again at 80 ° C for 1 min to improve adhesion by softening the CAB above its glass transition temperature.

To remove the CAB, the sample was first submerged in acetone for 15 min, followed by immersion in isopropyl alcohol (IPA) to eliminate acetone residues. Finally, a critical point dryer (CPD) was used: the sample, still submerged in IPA, was sealed in the CPD chamber, which was flooded with liquid CO 2 . The transition from supercritical CO 2 to gas prevented damage from surface tension. This process resulted in high-quality samples with minimal atomic-scale damage and reduced surface contamination.

Scanning transmission electron microscopy

A Nion UltraSTEM with beam energies between 60 keV and 200 keV was used with a nominal probe current of 20–60 pA, a 32 mrad semi-convergence angle and a collection angle of 80–200 mrad for HAADF imaging. Nanobeam electron diffraction (NBED) data were collected at 60 keV with a beam current of 50 pA and a total integration time of 5 s with a 0.3 mrad semi-convergence angle and a beam diameter of 50 nm. An Orca Hamamatsu CCD was used for collecting Ronchigrams. Image data were collected after defect patterning at 200 keV. The pressure in the objective area is near ultrahigh vacuum (<1 × 10 −9 mbar).

Custom beam trajectories were initially performed using a National Instruments multifunction I/O field-programmable gate array (FPGA, USB-7856) connected to the external beam motion input and HAADF output using National Instruments terminal block SCB-68At 54 . This platform operates on LabVIEW but commands are accessible by Python directly from the microscope user interface, Nion Swift.

A custom scan controller developed by AtomQ was used to add flexibility in applying different beam actions to generate user-selected defects within the same writing pattern, as shown in Fig. 1d, and to generate the large-area pattern in Fig. 5. The system integrates with Python and allows all beam control operations, such as pattern generation, parameter tuning and execution, to be performed directly from interactive notebooks.

In both cases, this allows changes from, for example, annular spiral beam motion parameters for the atomic lock-on algorithm (ALO; see next section), local positioning using the stay locked on algorithm (SLO), and repeated electron beam motions to steer Cr atoms.

The effectiveness of deterministic Cr displacement depends critically on the quality of the electron probe, sample tilt and defocus. Even small variations in probe aberrations such as astigmatism and spherical aberrations alter electron-wave propagation and channelling within individual atomic columns. As the displacement mechanism requires directing momentum and dose into a specific column with few-picometre lateral precision, subtle degradations in probe quality can lead to unintended scattering into neighbouring sites and reduce reproducibility. Therefore, maintaining a corrected and stable probe is essential for achieving the positioning accuracy and consistent Cr displacements demonstrated in this work.

Electron beam positioning

The level of control required for defect creation was enabled by developing multiple beam positioning and beam motion strategies 55 , 56 (Supplementary Information section 2 and Supplementary Figs. 2–6). The methods aim to expose the material only to electrons that actively contribute to atomic modification, minimizing unintentional structural changes. Target sites were defined as a sequence of 2D atomic vector steps. For the first target site, we carried out the ALO algorithm 53 , 55 , which performed a fast, sparse annular beam motion (100 ms duration, 1 nm radius) at a low dose rate (9.19 × 10 5 e − Å −2 s −1 ), achieving sub-20 pm lattice registration without the need for a dose-intensive overview image.

After precise positioning, the beam executed a pre-programmed beam motion (such as the ±b motion), then moved to the next target. For subsequent target positions, we used a second algorithm, SLO, which retrieved lattice coordinates with 25 pm accuracy in only 10 ms. This approach used 100 times fewer electrons than a single beam motion such as a 1 s long ±b beam motion. Between programmed targets, the electron beam was electrostatically blanked to eliminate unintended exposure during algorithmic computation. A detailed description of the technique used in this study will be provided in a separate, dedicated manuscript currently in preparation.

Atom steering

To control the displacement direction of Cr atoms within their atomic columns, we applied a directional beam motion of the electron beam at the target site 56 . The beam was moved repeatedly between two user-defined positions, a start and an end point. These positions were typically located on S/Br or Cr atomic columns, although arbitrary locations such as bond centres can also be used, as the beam trajectory is fully programmable. The repeated beam motion enabled real-time tracking of atomic configurations during electron beam motion by simultaneous signal acquisition on the HAADF detector. For all electron beam experiments presented here, a frequency of 0.5 kHz was used.

Atom streaking imaging

To probe and monitor Cr displacement dynamics, an imaging technique termed atom streaking was developed (Supplementary Information section 9). Analogous to streak camera imaging in energy- and time-resolved optical spectroscopy, this method enabled simultaneous spatial and temporal mapping of atomic configuration changes, with particular sensitivity to single-atom displacement dynamics.

A representative implementation of atom streaking involved applying a controlled, oscillatory electron beam motion (±b) across three atomic columns, as shown in Fig. 6. After achieving few-picometre precision in beam positioning using ALO or SLO, a HAADF detector signal was simultaneously acquired to monitor atomic rearrangements in real time. Each beam oscillation was integrated to form a single time step in the resulting image. Alternatively, with a sufficiently high signal-to-noise ratio, each individual beam motion was treated as a separate time step. The time resolution of this technique was directly determined by the beam oscillation frequency; in our typical setup at 0.5 kHz, one directional motion corresponded to 1 ms, whereas a complete oscillation (roundtrip) takes 2 ms. Higher frequencies allow for sub-millisecond time resolution.

This rapid, repeated, non-local beam control not only captured the displacement of atoms from their original positions but also captured their final locations. Owing to high detector efficiency, the technique provided unique insights into single-atom dynamics with few-picometre spatial precision and millisecond or even sub-millisecond temporal resolution, depending on the chosen beam motion frequency. In the example shown, beam motion was optimized to capture the transition of Cr atoms from their initial lattice site (δ0) to a nearby interstitial position (δ2). More generally, atom streaking is highly adaptable and was applied using all user-defined beam motion trajectories shown in Fig. 2.

Simulations

Image simulations

The intensity dependence of atomic columns was investigated using HAADF-STEM image simulations (Supplementary Information section 8 and Supplementary Fig. 24). Simulated 200 keV HAADF-STEM images were generated using the abTEM Python package 57 . Specifically, the simulations were performed using the PRISM algorithm 58 with a semi-convergence angle of 32 mrad, a collection angle range of 80–200 mrad, and an interpolation factor of 1. To match the experimental conditions, the simulated images were convolved with a Gaussian function with a full width at half maximum (FWHM) of 80 pm. Column intensities were extracted by averaging pixel values within a circular mask corresponding to the FWHM.

Electron beam channelling simulations

The effect of electron beam channelling (Supplementary Information section 3 and Supplementary Figs. 7–11) was investigated using the abTEM package 57 . The probe propagation through the sample was simulated using the multislice algorithm with a 200 keV electron beam and a semi-convergence angle of 32 mrad. The exit wave was computed and extracted at intervals of 0.99 Å, yielding a 3D representation of the beam intensity distribution throughout the sample. The code to replicate the simulations is available on GitHub ( ) with accompanying Python notebooks 59 .

Theoretical calculations

DFT calculations

For structure optimization and total energy calculations, we used the Vienna Ab Initio Simulation Package (VASP) 60 , 61 using projector augmented plane waves (PAW) 62 , 63 as the basis. To approximately address electronic correlations, we applied the rotationally invariant GGA+U approach with U = 2.5 eV and J = 0.4 eV (constrained random phase approximation values for single layer CrSBr in an effective dielectric environment of ε env ≈ 8; ref. 64 ) acting within the Cr d-orbitals. The PBE 65 form of the GGA exchange-correlation potentials was used. For the simulation of defects and their electronic structures, as well as the computation of single-particle (hopping t i j ) and many-body (bare Coulomb V i j k l ) matrix elements, we used N x × N y × 1 supercells with N x = 2, 3, 4, 5 and N y = 2, 3. We used Monkhorst–Pack k-grids of size K x × K y × 1, with K x and K y approximately inversely proportional to the corresponding in-plane lattice parameters of the supercell, and set a plane-wave cutoff energy of 500 eV. The migration energy of the Cr atom moving along the path between δ0 and δ2 sites was computed using a 4 × 3 × 2 supercell. All atomic positions were fully relaxed until the Hellmann–Feynman forces acting on each atom were reduced to 0.01 eV Å −1 or lower. Following DFT calculations, we constructed maximally localized molecular orbitals ψ using Wannier90 (ref. 66 ). On this basis, we calculated all single-particle (hopping) and many-body bare Coulomb matrix elements V i j k l . We estimated the partially screened density–density Coulomb matrix elements U i j as discussed in detail in Supplementary Information section 5.

Many-body perturbative bfQSboldsymbolGhatboldsymbolW approaches for computing electronic and excitonic spectra

The quasiparticle self-consistent GW approximation 37 , 38 is a self-consistent form of Hedin’s GW approximation. In contrast to conventional GW implementations, QSGW modifies the charge density and is determined by a variational principle 67 . Self-consistency in Σ ensures that feedback between the one-particle spectrum and the Coulomb interactions is accounted for 68 . This is particularly important when magnetic degrees of freedom are involved, as there is an additional coupling between spin and Σ. G, Σ and hatW are updated iteratively until all of them converge. Our results are thus parameter-free and have no starting point bias. Most of the discrepancies in estimating band gap and excitonic spectra in insulators originate from the omission of electron–hole interactions in the random phase approximation polarizability. By adding excitonic vertex (ladder diagrams) to the random phase approximation, we arrive at the QSGhatW approximation, in which electron–hole effects are taken into account. Generating hatW with excitonic ladder diagrams has several consequences. Most importantly, screening is enhanced and W is reduced. In the process of iterative self-consistency, this, in turn, reduces fundamental bandgaps and also valence bandwidths. The exact positions and wavefunctions of the defect states and the defect-bound excitonic states are sensitive to this mechanism of cross-coupling between spin and charge degrees of freedom in CrSBr. The theory has been shown to be successful in describing the electronic 69 , 70 and excitonic states in CrSBr 71 , 72 , 73 and its application to several weakly and strongly correlated magnetic insulators is given in refs. 36 , 74 . The importance of self-consistency in both QSGW and QSGhatW for different materials have been explored 75 .

For CrSBr, we simulated a 3 × 4 × 1 supercell containing 72 atoms within QSGhatW with individual layers containing ferromagnetically polarized spins pointing along the b-axis, while the magnetic interlayer coupling between adjacent layers was kept antiferromagnetic. The single-particle calculations (LDA and energy band calculations with the static quasiparticle self-energy obtained within QSGW and QSGhatWΣ(k)) are performed on a 4 × 4 × 8 k-mesh while the relatively smooth dynamical self-energy Σ(ω) is constructed using a 2 × 2 × 4 k-mesh. The QSGW and QSGhatW cycles are iterated until the RMS change in the static part of the quasiparticle self-energy Σ(0) reaches 10 −5 Ry. The two-particle Hamiltonian that is solved self-consistently to compute both the Σ and the excitonic eigenvalues and eigenfunctions, contained 108 valence bands and 38 conduction bands. For all defect simulations, the 1.3 eV excitonic state from the pristine compound still appears along with several defect-bound excitonic states. For example, for the V Cr +Cr int delta 2 complex, the deepest lying one-particle defect state appears at 0.7 eV, whereas the defect-bound excitonic state (with largest binding energy) appears at 0.47 eV, suggesting strong electron–hole binding mediated by the defect state.

Data post-processing

No additional post-processing steps of the HAADF-STEM image data were used to generate the figures, which are shown as-acquired and unfiltered.