Yoonseok Baek

Hilton B De Aguiar

Sylvain Gigan

Phase conjugation with spatially incoherent light in complex media

Shaping light deep inside complex media, such as biological tissue, is critical to many research fields. Although the coherent control of scattered light via wavefront shaping has made significant advances in addressing this challenge, controlling light over extended or multiple targets without physical access inside a medium remains elusive. Here we present a phase conjugation method for spatially incoherent light, which enables the non-invasive light control based on incoherent emission from multiple target positions. Our method characterizes the scattering responses of hidden sources by retrieving mutually incoherent scattered fields from speckle patterns. By time-reversing scattered fluorescence with digital phase conjugation, we experimentally demonstrate focusing of light on individual and multiple targets. We also demonstrate maximum energy delivery to an extended target through a scattering medium by exploiting transmission eigenchannels. This paves the way to control light propagation in complex media using incoherent contrasts mechanisms.

I. Main

Delivering optical energy and transmitting information through complex media remains an important challenge in many fields of studies, including optical manipulation [1], deep-tissue imaging [2,3] and optogenetics [4,5]. In recent years, it has been shown that the coherent control of scattered light can manipulate spatial, spectral and temporal distributions of light in scattering media [6][7][8][9][START_REF] Cao | Shaping the propagation of light in complex media[END_REF]. However, such capabilities are greatly limited without physical access inside a medium because the scattering response to target position is difficult to characterize. As a result, non-invasive light control over extended or multiple targets remains elusive despite being crucial for real-world applications. Optimizing incident wavefront based on a feedback signal [START_REF] Katz | Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers[END_REF][START_REF] Daniel | Light focusing through scattering media via linear fluorescence variance maximization, and its application for fluorescence imaging[END_REF][START_REF] Boniface | Noninvasive light focusing in scattering media using speckle variance optimization[END_REF][START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF][START_REF] Li | Non-invasive optical focusing inside strongly scattering media with linear fluorescence[END_REF][START_REF] Rauer | Scattering correcting wavefront shaping for three-photon microscopy[END_REF][START_REF] Thompson | Wavefront shaping enhanced raman scattering in a turbid medium[END_REF][START_REF] Tian | Noninvasive chemically selective energy delivery and focusing inside a scattering medium guided by raman scattering[END_REF] is mostly limited to focusing on a single isolated target, and even then it has limitations that require numerous changes of the wavefront. While time-reversal or phase conjugation techniques [START_REF] Yaqoob | Optical phase conjugation for turbidity suppression in biological samples[END_REF][START_REF] Cui | Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation[END_REF][START_REF] Hsieh | Digital phase conjugation of second harmonic radiation emitted by nanoparticles in turbid media[END_REF][START_REF] Xu | Timereversed ultrasonically encoded optical focusing into scattering media[END_REF][START_REF] Ivo | Digital optical phase conjugation of fluorescence in turbid tissue[END_REF][START_REF] Judkewitz | Speckle-scale focusing in the diffusive regime with time reversal of variance-encoded light (trove)[END_REF][START_REF] Ma | Time-reversed adapted-perturbation (trap) optical focusing onto dynamic objects inside scattering media[END_REF][START_REF] Haojiang Zhou | Focusing on moving targets through scattering samples[END_REF][START_REF] Haowen Ruan | Optical focusing inside scattering media with time-reversed ultrasound microbubble encoded light[END_REF][START_REF] Ruan | Focusing light inside scattering media with magnetic-particle-guided wavefront shaping[END_REF][START_REF] Yang | Focusing light inside live tissue using reversibly switchable bacterial phytochrome as a genetically encoded photochromic guide star[END_REF][START_REF] Aizik | Fluorescent wavefront shaping using incoherent iterative phase conjugation[END_REF] allow for effective light delivery to an optical or virtual source, they cannot individually control light on multiple targets.

Here we address these challenges by extending phase conjugation to spatially incoherent light. Our approach utilizes incoherent emission from multiple targets. We first characterize scattering responses of these hidden sources by retrieving mutually incoherent fields from spatially modulated speckle patterns. The retrieved fields are related to the field transmission matrix [START_REF] Sébastien | Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media[END_REF], and their phase conjugation enables light control over the desired positions. We demonstrate this experimentally by focusing light on individual and multiple fluorescent targets through a scattering medium. Finally, we show that transmission eigenchannels can be identified by decomposing the incoherent fields and demonstrate maximum energy delivery to a hidden extended target. * yoonseok.baek@lkb.ens.fr

II. Results

A schematic of the experiment is illustrated in Fig. 1. We consider a scenario where multiple fluorescent targets are hidden by a scattering medium and act as guidestars. Fluorescence emitted by these guidestars is scattered, resulting in an incoherent addition of speckle patterns on the camera. Under this condition, we aim to deliver light back to each of the guidestars by time-reversing the scattered fluorescence. Our approach consists of retrieving multiple incoherent fields that compose the fluorescence, and using them to generate phase-conjugated beams. To this end, we introduce wavefront modulation of the scattered fluorescence with a spatial light modulator (SLM). This modulation induces changes in the measured incoherent speckle patterns, providing information to retrieve the scattered fields that will later be used for phase conjugation. To explain the retrieval process in detail, we introduce a partial field transmission matrix T whose input and output are fields at the SLM plane and guidestar positions, respectively, E guidestar = T E SLM . According to the time-reversal symmetry, the scattering of fluorescence emitted by N guidestars at the same wavelength can be expressed by the rows of the transmission matrix, T = [t 1 , . . . , t N ] ⊤ , where t n represents the scattered field that corresponds to an individual guidestar. For the sake of simplicity, we will assume that the brightness of each guidestar is the same (see Supplementary Information Section 6 for general cases with different brightness). Then, the incoherent speckle pattern measured by the camera is expressed as I 0 = diag (T F ) † T F , where F represents the discrete Fourier transform between the SLM and camera plane (see Supplementary Information Section 1 for the derivation). If we apply M different wavefront modulations using the SLM, and express the mth modulation as a diagonal matrix S m , the incoherent speckle pattern after the modulation is

I m = diag (T S m F ) † T S m F . (1)
To find a set of scattered fields that satisfy Eq. used an iterative approach that minimizes the error between the measured and predicted I m . We began by making an initial guess of t n and applying the wavefront modulation numerically according to Eq. 1. Then, following the mixed-state reconstruction [START_REF] Thibault | Reconstructing state mixtures from diffraction measurements[END_REF], we corrected the amplitude parts of the modulated speckle fields at the camera plane, using the ratio of the measured to the predicted I m . Next, we compensated for the wavefront modulation and updated the guess of t n . These steps were repeated for all the wavefront modulations, resulting in a maximum likelihood estimation of T (see Supplementary Information Section 2 for more information). We note that there is inherent ambiguity in determining T , as T and its unitary transformation are indistinguishable based on intensity. This can be confirmed by replacing T with U T in Eq. 1. For this reason, the scattered fields are retrieved as a mixture of t n :

H = U T , ( 2
)
where H is a set of retrieved fields h n , H = [h 1 , ..., h N ] ⊤ , and U is an arbitrary unitary matrix. We further note that Eq. 2 arises from the nature of mutual incoherence, rather than from the reconstruction method. Despite this ambiguity, the retrieved fields offer unique capabilities for phase conjugation, as we show below.

Incoherent phase conjugation

The time-reversal of the entire scattered fluorescence will regenerate light at hidden sources, creating foci on the entire targets. One way to accomplish this is through the ensemble average of the phase-conjugated scattered fields t * n . Alternatively, we chose to use the ensemble average of h * n because it gives the identical phase conjugation result (see Supplementary Information Section 3). In both cases, the ensemble average results in the incoherent sum of the phase conjugation of the scattered fields, which can be realized by shaping either coherent or incoherent light. We will refer to this approach as "Incoherent phase conjugation" for simplicity, regardless of the coherence of light used to generate the phase-conjugated fields (see Discussion for more information).

To demonstrate this incoherent phase conjugation, we introduced several 1 µm fluorescent beads as guidestars behind the scattering medium. We retrieved the multiple scattered fields according to the number of the guidestars. In our experiments, each 1 µm bead was considered as an individual guidestar as the speckle grain size at the target plane was ∼0.9 µm. To implement the incoherent phase conjugation, we used the SLM and a collimated laser beam to generate N different phase-conjugated fields of h * n over time and measured their time-averaged response for phase (see Methods and Discussion). To evaluate the performance of the phase conjugation, we first conducted an experiment with a single fluorescent bead (Fig. 2a). By phase-conjugating the scattered field, we observed a strong focus on the bead (Fig. 2b). This is in clear contrast with the random speckle generated by a beam with a random wavefront (Fig. 2c). The enhancement factor, defined as the ratio between the optimized focus intensity and mean background intensity, was ∼4,400. We then placed multiple fluorescent beads (Fig. 2d-f) behind the scattering medium. By incoherently phase conjugating the scattered fields, we successfully generated foci at every guidestar positions (Fig. 2g-i). Despite the minimal spectral memory effect [START_REF] Vesga | Focusing large spectral bandwidths through scattering media[END_REF], we were also able to excite the bead through the scattering medium by shaping the excitation beam using the retrieved field at the emission wavelength(Fig. S1).

Selective focusing on individual targets

In order to selectively focus on individual targets, it is necessary to demix the individual fields t n from their mixture h n . We note that t n is not strictly orthogonal, and thus the orthogonalization of h n can not be a solution. Our solution was to directly invert Eq. 2 by finding U . To this end, we utilized the memory effect [START_REF] Freund | Memory effects in propagation of optical waves through disordered media[END_REF][START_REF] Feng | Correlations and fluctuations of coherent wave transmission through disordered media[END_REF], where neighboring guidestars generate correlated speckle patterns. Specifically, we iteratively applied a random unitary transformation to the retrieved fields H, such that the correlation between the transformed speckle patterns is maximized (see Supplementary Information Sec- tion 5).

Figure 3 shows the experimental result with 5 fluorescent beads. When the scattered fields h n are directly used for phase conjugation, each phase conjugation generated foci on several guidestars with different intensities (Fig. 3b). The demixed fields, on the other hand, generated a focus on a single guidestar, showing that t n is successfully recovered (Fig. 3c). The memory effect range in this experiment, defined as the full width at half maximum of speckle cross-correlation, was 5 µm, which is much smaller than the spatial extent of the guidestars. This result shows that the selective focusing via demixing is possible as long as a pair of guidestars lies within the memory effect range.

Targeted energy delivery

Maximum energy delivery through scattering media requires an eigenchannel of T [START_REF] Ivo | Universal optimal transmission of light through disordered materials[END_REF][START_REF] Kim | Maximal energy transport through disordered media with the implementation of transmission eigenchannels[END_REF]. The transmission eigenchannels correspond to the singular vectors of T , and the first singular vector with the largest singular value, delivers the maximum energy to the target. Although the direct access to T is not always possible, the transmission eigenchannels of T can be found using H. This is because the eigenchannels of T and H are identical because H † H = T † T according to Eq. 2. Thus, we can deliver the maximum energy to extended targets using the first singular vector of H.

To demonstrate the targeted energy delivery, we placed a 5 µm fluorescent ink droplet behind the scattering medium (Fig. 4a). Based on the size of the target, we estimated the number of incoherent fields and retrieved 23 scattered fields. We then performed the singular value decomposition of the retrieved fields H. Finally, we injected fields that corresponds to the singular vectors v n and observed the energy delivered to the target. When a random wavefront is injected to the scattering medium, a speckle pattern is generated at the target plane (Fig. 4b). In contrast, the singular vectors produce intensity distributions highly concentrated on the target (Fig. 4c). By summing the results of all the singular vectors, we confirmed that the energy is delivered only to the target area (Fig. 4d). The first singular vector v 1 shows a 174-times increase in the energy on the target, compared to random wavefronts. The enhancement decays with the singular vector index [Fig. 4(e)]. We observed that the values are not perfectly sorted in a descending order, which we believe is due to the numerical error in the retrieved fields and to the use of the phase-only SLM.

III. Discussion and conclusion

We have presented an approach to control light in scattering media without the physical access to a target plane by extending phase conjugation to spatially incoherent light. We have demonstrated focusing and maximum energy delivery to extended fluorescent targets through a scattering medium. Our method, like other digital phase conjugation or wavefront shaping techniques, can provide reliable focusing results with a high contrast provided a sufficiently large number of incident modes are controlled [7]. Nonetheless, what distinguishes our approach is its ability to tackle mutually incoherent fields of scattered fluorescence. Another important aspect is that it does not require the precise alignment between the camera and SLM [START_REF] Timothy R Hillman | Digital optical phase conjugation for delivering two-dimensional images through turbid media[END_REF], since the scattered fields are retrieved at the SLM plane. Its principle of characterizing the scattering response is entirely passive, as it does not alter the emission of guidestars, as opposed to to techniques that modulates the excitation wavefront (e.g. [START_REF] Boniface | Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix[END_REF]). It only requires that the emission from the sources be constant on average over the measurement of I m and thus unaffected by the setup geometry.

To further demonstrate the capability of the proposed method under realistic conditions, we conducted an additional experiment. Specifically, we employed a volumetric scattering medium composed of parafilm layers, and we increased the number of guidestars. More information on the experiment can be found in Supplementary Information Section 4. The results showed that the proposed method was able to produce high-contrast foci on 13 and 22 target positions (Fig. S2), demonstrating the potential of our method for practical applications.

For the proposed method, it is important to estimate the number of mutually incoherent waves, N . This is because the underestimation of N results in imperfect reconstruction of the scattered fields. We note that the overestimation is allowed because it results in redundant reconstruction (see Fig. S3). Nevertheless, it is recommended to use the smallest possible value of N for the minimal measurements and computation time. There are different methods to estimate N . The contrast of fluorescence speckle is an useful indicator for N , as it decreases as √ N [START_REF] Joseph W Goodman | Speckle phenomena in optics: theory and applications[END_REF]. It is also possible to find N by analyzing the error in the mixed-state phase retrieval or the singular value distribution [START_REF] Michel Popoff | Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern anal-ysis[END_REF] for different values of N . In our demonstrations, we did not consider the spectral degrees of freedom because narrow spectral responses were measured by using interference filters. If the detection bandwidth is greater than the spectral memory effect range, different spectral components should be considered in estimating N . This will allow the proposed technique to handle situations where hidden sources emit at different wavelengths.

Another important consideration is the number of modulation M required for the field retrieval, which scales linearly with N . In experiments, reliable phase conjugation results were obtained when M ≥ 6N (see Supplementary Information Section 2). This linearity can be attributed to the multiplexed information in the intensity of multiple incoherent fields. We emphasize that M scales with the number of incoherent waves N not with the number of controlled input modes of the SLM. Recent advances in phase retrieval [START_REF] Emmanuel | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Emmanuel | Phase retrieval via matrix completion[END_REF] show that few measurements are sufficient in retrieving a coherent field (N = 1). In this regard, we believe that in principle even fewer M may be used for our method.

In our proof-of-principle experiments, we used the simple algorithms for the retrieval and demixing of the scattered fields. The performance of the algorithms can be enhanced by incorporating constraints, convex optimization [START_REF] Emmanuel | Phase retrieval via matrix completion[END_REF] or the generalized memory effect [START_REF] Osnabrugge | Generalized optical memory effect[END_REF]. The total measurement time can be reduced by designing a setup with minimal energy loss, and by using a sensitive detector, such as an EMCCD.

The current implementation of incoherent phase conjugation is based on the alternation of multiple coherent fields. This method can accurately generate the incoherent response of the time-reversed fluorescence, although it requires some integration time [START_REF] De Santis | Synthesis of partially coherent fields[END_REF]. Importantly, the required integration time can be reduced to few microseconds by employing fast spatial light modulators [START_REF] Rodenburg | Experimental generation of an optical field with arbitrary spatial coherence properties[END_REF][START_REF] Tzang | Wavefront shaping in complex media with a 350 khz modulator via a 1d-to-2d transform[END_REF]. Alternatively, the incoherent phase conjugation can be implemented by directly shaping an incoherent source [START_REF] Barré | Holographic beam shaping of partially coherent light[END_REF]. In this case, the incoherent response will be obtained at a timescale greater than the coherence time of the source.

In conclusion, our method enables versatile light control over extended or multiple targets using incoherent contrast mechanisms. The concept can be applied to different incoherent emissions, such as spontaneous Raman scattering [START_REF] Thompson | Wavefront shaping enhanced raman scattering in a turbid medium[END_REF][START_REF] Tian | Noninvasive chemically selective energy delivery and focusing inside a scattering medium guided by raman scattering[END_REF], and a wide range of photoluminescence [START_REF] Kim | Quantum reference beacon-guided superresolution optical focusing in complex media[END_REF]. Furthermore, it enables the passive characterization of a transmission matrix, opening up the possibility to generalized light control using transmission-matrix-based operators [START_REF] Cao | Shaping the propagation of light in complex media[END_REF]. We envision that the proposed approach will enable targeted light delivery through thick biological tissue, facilitating biomedical applications, such as optogenetic stimulation and phototherapy.

IV. Methods

Experimental setup

The experimental setup is shown in Fig. S4. A laser diode (λ = 488 nm, LP488-SF20G, Thorlabs) was used to excite fluorescence. The excitation beam was delivered to guidestars by a lens (L1, f = 100 mm) and an objective lens (Plan N 20× 0.4, Olympus). To monitor the guidestars and phase conjugation, a dichroic mirror (DMLP490R, Thorlabs), a lens (L2, f = 200 mm), a bandpass filter (FL532-10, Thorlabs), and a camera (acA5472-17um, Basler) were placed on the side without a scattering medium. The guidestars were fluorescent ink mixed with UV glue (NOA 68, Norland), and fluorescent beads (F8823, Invitrogen) immersed in glycerol. A scattering medium was a 220-grit ground glass diffuser, placed approximately 170 µm away from the guidestars. On the detection side, the scattered fluorescence was collected by an objective lens (MPlan N 50× 0.75, Olympus) and two lenses (L3, f = 75 mm; L4, f = 150 mm). An SLM (X10468-04, Hamamtsu) and a linear polarizer were used to modulate the fluorescence. The modulated fluorescence is Fourier transformed by lenses (L5, f = 100 mm; L6, f = 200 mm; L7, f = 250 mm) and then captured by an sCMOS camera (PCO.edge 5.5, PCO) with bandpass filters (BP2, FL532-3 and FBH520-40, Thorlabs). For the fluorescent beads we used, these filters transmit 5.35% of the total fluorescence signal. An iris was placed between L6 and L7 to adjust the speckle grain size at the camera. For phase conjugation, a laser (λ = 532 nm, Compass 215M-50, Coherent) was collimated using a 5 µm pinhole and a lens (L8, f = 6 mm). The collimated laser was then shaped by the SLM to generate a phaseconjugated beam, which propagated back to the medium. A flip mirror was used to switch between the fluorescence detection and phase conjugation.

Phase-conjugated beam generation

A phase-conjugated beam is generated using the collimated laser beam and the SLM. The collimated beam is shaped to the phase conjugate of a given scattered field E scattered by displaying a phase pattern that corresponds to -arg (E scattered ) on the SLM. The resultant phaseconjugated beam propagates back through the scattering medium, retracing the scattering paths of fluorescence.

where h nk is the kth element of h n . In our experiments, we stopped the algorithm when the MSE was less than 10 -5.6 . See Fig. S6 for a change of the MSE over iterations. We confirmed numerically that the algorithm retrieves a correct set of fields for M ≥ 4N in the absence of measurement noise. The size of macro-pixel had almost no effect on the reconstruction, except when its size is comparable to the SLM. We observed that the minimum value of M required for correct reconstruction increases depending on the noise level. For the experimental results shown in the main text, we used 6-8N modulations.

Incoherent phase conjugation

When all the mutually incoherent components of fluorescence are time-reversed, the intensity at the nth guidestar can be expressed by multiplying the nth row of the transmission matrix, t ⊤ n , and phase-conjugated field t * m :

m t ⊤ n t * m 2 = t ⊤ n T † T t * n , (S5)
where * denotes the complex conjugate. For a medium with negligible reflection and absorption, T † T ≈ 1. Thus Eq. S5 is simplified to t ⊤ n t * n , which remains more or less constant regardless of n. As a result, the incoherent phase conjugation of t n generates foci on the entire guidestars with roughly the same intensity. It should be noted that T is not unitary in practice, since we are not able to measure the entire scattered light. Nevertheless, if the loss of information due to absorption or measurement is minimal, the phase conjugation can approximate the case of unitary T for complex media, producing high-contrast foci on the entire guidestars [4,5]. Similarly, the intensity for the incoherent phase conjugation of h n is expressed as,

m t ⊤ n h * m 2 = t ⊤ n H † Ht * n . (S6)
Since H † H = T † T according to Eq. 2, Eq. S6 is identical to Eq. S5. Thus, the incoherent phase conjugation of h n and t n are identical.

Additional experiment with volumetric scattering medium

In addition to the incoherent phase conjugation described in Fig. 2, we conducted an additional experiment with a volumetric scattering medium. In the experiment, we placed 1 µm fluorescent beads directly behind three layers of parafilm, with a total thickness of approximately 360 µm corresponding to half of the transport mean free path [6]. We replaced the objective lens with a long-working-distance objective lens (MPlanFL N 100× 0.9, Olympus) and used a 10-nm bandpass filter (FLH05532-10, Thorlabs) which transmits 14.3% of the fluorescence signal. The rest of the setup was the same. The experiment was conducted in the same manner as the experiment for Fig. 2.

The results showed that that the proposed method was able to produce high-contrast foci on 13 and 22 guidestars (Fig. S2). In Fig. S2d, we observed that foci were not generated at a few target positions located far from the center. This is due to the limited field-of-view (FoV) in which the SLM can control light. By evaluating the speckle pattern generated by a back-propagating beam with a random wavefront, we confirmed that the size of the controllable FoV was approximately 60 µm (Fig. S2e). This is smaller than the field of view of Fig. S2 (75.6 µm). We note that the FoV is determined by the scattering media and the experimental setup, such as the focal length of an objective lens and the SLM pixel size.

Demixing incoherent fields

The scattered fields t n can be recovered by reversing Eq. 2. To this end, we developed an algorithm inspired by the simulated annealing [7] that exploits the memory effect between the speckle patterns at the camera plane. The algorithm starts by applying a random unitary transformation to the retrieved fields H and expressing the result at the camera plane:

H = U j fj Û HF , ( S7
)
where j is the iteration index, U j is a random unitary matrix, f j is an exponent that decays to 0 over the iteration, and Û is the current best estimate of U -1 . The decaying behavior of f j makes U j induce smaller changes as the iteration progresses: lim j→∞ U j fj = 1. Then the correlation metric C is calculated at every iteration,

C = N n=1 max |h n | 2 ⋆ |h (n mod N )+1 | 2 2 , ( S8
)
where h n (k) is the unitary-transformed speckle field at the camera plane, corresponding to the nth row of H, and ⋆ is the cross-correlation. If the value of C is greater than the previous values, Û is updated to U j fj Û . If not, Û is unchanged. At the end of the iteration, the algorithm results in Û ≈ U -1 , and the scattered fields are recovered by inverting Eq. 2: T = Û H. We initialized the iteration with Û = 1 and f j = 100/(j + 1), and the maximum iteration number of 10 4 . We note that the experimental noise in H can lead to an incorrect demixing. To mitigate this issue, we calculate C only using speckles brighter than the average, |h n | 2 > µ, where µ is the mean value of every |h n | 2 .

Guidestars of varying brightness

If the brightness of the guidestars is different, we can express the scattered field corresponding to nth guidestar as b n t n , where b n is a constant that is the square-root of the energy emitted by the guidestar and t n is a normalized field, such that t † n t n = 1. Then we can rewrite the expressions in the main text by replacing T with BT , where B = diag(b 1 , ..., b n ). For example, the unmodulated speckle image becomes I 0 = diag (BT F ) † BT F . Similarly, Eq. 1 is rewritten as

I m = diag (BT S m F ) † BT S m F , (S9)
and Eq. 2 as

H = U BT . (S10)
We note that this generalization is related to the interpretation of the measured image and retrieved fields. Therefore, the mixed-state reconstruction and the demixing process are unaffected by the brightness of the guidestars. For incoherent phase conjugation (Supplementary Section 3), foci will be generated on the entire guidestar with the intensities that correspond to the brightness of the guidestars. Regarding maximum energy delivery, the eigenchannel of H will preferentially couple light into the brighter guidestar positions because

H † H = T † B 2 T , where B 2 = diag(b 2 1 , ..., b 2 n ).

Maximum number of guidestars

In principle, the mixed-state algorithm can find solutions regardless of the number of guidestars N , provided that there are enough wavefront modulations. However, in practice, the number of guidestars could be limited because the dynamic range of a camera causes quantization of images. For the phase retrieval of coherent light, it is known that 3-bit representation of speckle can produce reliable phase retrieval results [8]. Considering that the mixed-state reconstruction can be thought of as the information multiplexing of mutually incoherent signals, we would expect the following condition to be satisfied,

N × 2 3 ≤ 2 (Camera bit depth) . (S11)
For example, with a 16-bit resolution, the maximum theoretical N will be 8192, although in practice it would be significantly lower. Finally, we note that even a larger number could be allowed by incorporating the quantization error into the reconstruction algorithm [9].

Figure 1 .

 1 Figure 1. Schematic of phase conjugation with incoherent fluorescence: Multiple fluorescent guidestars are hidden behind a scattering medium. The scattered fluorescence, composed of mutually incoherent speckle fields, is modulated by an SLM. The modulated fluorescence is Fourier transformed by a lens and results in the incoherently added speckle patterns on the camera. Then mixed-state phase retrieval recovers a set of scattered fields, whose phase conjugation enables targeted light control.

Figure 2 .

 2 Figure 2. Incoherent phase conjugation for multiple targets: a, The fluorescence image of an 1 µm bead taken from the side without a scattering medium. b, Intensity at the target plane generated by phase conjugation based on the scattered fluorescence field. (Inset) The phase of the scattered field shown in the HSV colormap. The central highlighted part is used for the phase conjugation. c, Intensity at the target plane generated by a random wavefront. d-f, The fluorescence images of multiple beads hidden behind the scattering medium. g-i, Intensity at the target plane with incoherent phase conjugation of scattered fields.

Figure 3 .Figure 4 .

 34 Figure 3. Selective focusing on individual targets: a, The fluorescence image of a target comprised of 1 µm beads. b,c, Intensity at the target plane after the phase conjugation of individual scattered fields before (b) and after (c) the demixing process. (Insets) The phase of the scattered fields used for phase conjugation.

Figure S1 .

 S1 Figure S1. Phase conjugation using an excitation beam: a Fluorescence image of a bead when an excitation beam (λ = 475 µm) is shaped using the phase conjugation pattern of scattered fluorescence (λ = 532 µm). b Fluorescence image of the same bead when the excitation beam is shaped by a random phase pattern. c Intensity profiles of (a) and (b) along the horizontal line that crosses the center.

Figure S2 .Figure S3 .Figure S4 . 1 Figure S5 .

 S2S3S41S5 Figure S2. Incoherent phase conjugation in volumetric scattering medium: a-b Fluorescence image of beads hidden behind 3 layers of parafilm. c-d Results of the incoherent phase conjugation. e Intensity at the target plane with a random wavefront.

Mixed-state phase retrieval Wavefront synthesis

	Scattering medium	SLM	Lens	Camera images
	Fluorescent guidestars			1	2	...	M
	1	2	M ...
				Scattered fields
						1	2	.. .	N

1, we arXiv:2301.13140v2 [physics.optics] 29 May 2023

Acknowledgement

This research was supported by H2020 Future and Emerging Technologies (863203), European Research Council (724473), and Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2022R1A6A3A03072108).

Supplementary Information: Phase conjugation with spatially incoherent light in complex media 1 Matrix representation of I 0 and I m

The intensity speckle measured at the camera plane can be written as

where t n (r) and tn (k) is the nth scattered field at the SLM plane and camera plane, respectively, and F is the Fourier transform. To express Eq. S1 in a matrix form, we utilize T to express t n (r) and tn (k). First, we can express t n (r) as the rows of T , such that T ij = t i (r j ). Similarly, tn (k) can be expressed as the rows of T F , where F is a Fourier transform matrix, such that (T F ) ij = ti (k j ). If we consider a matrix (T F ) † T F , its elements can be expressed as

Finally, we can rewrite Eq. S1 as I 0 = diag (T F ) † T F by expressing the main diagonal of a matrix using the function diag. Equation 1 can be obtained by following a similar process as described above, but by replacing T with T S m .

2 Mixed-state phase retrieval

Our phase-retrieval algorithm is designed to find N mutually incoherent fields at the SLM plane from M modulated intensity images. We initialize the algorithm by letting the scattered fields h n as N complex Gaussian random fields. Then we apply the wavefront modulation S m to the scattered fields h n . We note that S m is a diagonal matrix whose diagonal elements correspond to the field modulation given by the SLM. In the experiments, SLM was divided into macro-pixels (composed of 40 × 40 pixels) having random phase values (0 or π). After the wavefront modulation S m , the field at the camera plane is expressed as h(m) n = (S m F ) ⊤ h n . Then we conduct Fourier magnitude projection using an auxiliary function ψ:

where k is a coordinate in the spatial frequency domain, h(m) n (k) is the modulated field at the camera plane, corresponding to h(m) n , and γ is a constant parameter. Then the fields are updated by compensating the wavefront modulation:

where

n (k). The update through Eq. S2-S3 is continued for the entire M modulations. The whole process is repeated for several times to obtain consistent h n (see Fig. S5). With γ = 1/2, this method can be interpreted as the maximum likelihood reconstruction and retrieves the incoherent fields [1]. However, we observed that few initial iterations with γ = 1 accelerates the convergence greatly, in both numerical simulations and experiments. This observation agrees with the effect of the nonlinear modulus constraint that was reported in [2,3]. It should be noted when γ is greater than 1/2, the brighter area of I m carry more weight during the iteration. While it has been reported that the larger values of γ lead to faster convergence [2], we expect that the optimal value of γ for the initial iterations would depend on the specific application. In our experiments, we used γ = 1 for the first 20 iterations and γ = 1/2 for the rest of the iterations. We stopped the algorithm when there was no significant change in h n . Specifically, we evaluated the convergence by using the mean squared error (MSE) between the current and previous estimates of h n :


Arjun Ashoka

Adam J Clancy

Naitik A Panjwani #

Nicholas J M Popiel

Alex Eaton

Thomas G Parton

Loren Picco

Sascha Feldmann

Rebecca R C Shutt

Remington Carey

Eva S Y Aw

Thomas J Macdonald

Marion E Severijnen

Sandra Kleuskens

Hilton Barbosa De Aguiar

Richard H Friend

Jan Behrends

Peter C M Christianen

Christopher A Howard

Akshay Rao

Raj Pandya

Room Temperature Optically and Magnetically Active Edges in Phosphorene Nanoribbons 1

 . Here, we report the magnetic and semiconducting properties of phosphorene nanoribbons (PNRs). Static (SQUID) and dynamic (EPR) magnetization probes demonstrate that at room temperature, films of PNRs exhibit macroscopic magnetic properties, arising from their edge, with internal fields of ~ 250 to 800 mT. In solution, a giant magnetic anisotropy enables the alignment of PNRs at modest sub-1T fields. By leveraging this alignment effect, we discover that upon photoexcitation, energy is rapidly funneled to a dark-exciton state that is localized to the magnetic edge and coupled to a symmetry-forbidden edge phonon mode. Our results establish PNRs as a unique candidate system for studying the interplay of magnetism and semiconducting ground states at room temperature and provide a stepping-stone towards using low-dimensional nanomaterials in quantum electronics.

Over the last 40 years, strategies to couple different magnetic phenomena and the electronic ground states of semiconductors have mostly focussed on doping traditional semiconductors with certain transition metals to form so-called dilute magnetic semiconductors (DMSs) [START_REF] Náfrádi | Optically switched magnetism in photovoltaic perovskite CH3NH3(Mn:Pb)I3[END_REF][START_REF] Story | Carrier-concentration-induced ferromagnetism in PbSnMnTe[END_REF][START_REF] Haury | Observation of a Ferromagnetic Transition Induced by Two-Dimensional Hole Gas in Modulation-Doped CdMnTe Quantum Wells[END_REF][START_REF] Wunderlich | Coulomb blockade anisotropic magnetoresistance effect in a (Ga,Mn)As single-electron transistor[END_REF][START_REF] Moreno | Correlation of structure and magnetism in GaAs with embedded Mn(Ga)As magnetic nanoclusters[END_REF] . A more recent approach has been to exploit the unique spin configurations that can be realised in two-dimensional (2D) monolayer materials -some of which, such as CrSBr and NiPS3, have been also shown to host (at low temperature) semiconducting magnetic phases [START_REF] Bae | Exciton-coupled coherent magnons in a 2D semiconductor[END_REF][START_REF] Lee | Magnetic Order and Symmetry in the 2D Semiconductor CrSBr[END_REF][START_REF] Wang | Spin-induced linear polarization of photoluminescence in antiferromagnetic van der Waals crystals[END_REF][START_REF] Kang | Coherent many-body exciton in van der Waals antiferromagnet NiPS 3[END_REF] . Slicing two-dimensional monolayers into thin strips -nanoribbons -opens the possibility for an even wider range of magnetic spin-arrangements across and along the ribbon edges, whilst retaining the favourable electronic and mechanical properties of 2D systems. The most intensely studied nanoribbons thus far have been those based on the flat, honeycomb lattice of graphene [START_REF] Magda | Room-temperature magnetic order on zigzag edges of narrow graphene nanoribbons[END_REF][START_REF] Zhao | Topological Phases in Graphene Nanoribbons Tuned by Electric Fields[END_REF] , where at low temperatures topologically-engineered bands [START_REF] Rizzo | Inducing metallicity in graphene nanoribbons via zero-mode superlattices[END_REF][START_REF] Rizzo | Topological band engineering of graphene nanoribbons[END_REF][START_REF] Gröning | Engineering of robust topological quantum phases in graphene nanoribbons[END_REF] and a range of spin ordered electronic states have been observed [START_REF] Slota | Magnetic edge states and coherent manipulation of graphene nanoribbons[END_REF][START_REF] Blackwell | Spin splitting of dopant edge state in magnetic zigzag graphene nanoribbons[END_REF] . However, the laborious chemical synthesis [START_REF] Cai | Atomically precise bottom-up fabrication of graphene nanoribbons[END_REF][START_REF] Ruffieux | On-surface synthesis of graphene nanoribbons with zigzag edge topology[END_REF] combined with challenging edge-functionalisation [START_REF] Slota | Magnetic edge states and coherent manipulation of graphene nanoribbons[END_REF][START_REF] Blackwell | Spin splitting of dopant edge state in magnetic zigzag graphene nanoribbons[END_REF] -akin to magnetic doping in DMSsnecessary to achieve the aforementioned phenomena, represents a significant challenge for both practical use and exploring novel physics in GNRs. These problems are further compounded by the short length (sub-100 nm) that GNRs can be limited to and the fact that spintronic properties often emerge only at low-temperatures. Consequently, there is a need to explore nanoribbon systems which can be produced at scale and exhibit room-temperature intrinsic spintronic and electronic properties. Among the mooted systems, phosphorene nanoribbons (PNRs) -the black phosphorus analogue to GNRs -have uniquely been proposed. Although PNRs are expected to host many superlative characteristics, the predicted presence of excitonic effects and room temperature ferromagnetism [START_REF] Wu | Spin density waves predicted in zigzag puckered phosphorene, arsenene and antimonene nanoribbons[END_REF][START_REF] Sisakht | Strain-induced topological phase transition in phosphorene and in phosphorene nanoribbons[END_REF][START_REF] Nourbakhsh | Excitons and optical spectra of phosphorene nanoribbons[END_REF][START_REF] Wu | Electronic and transport properties of phosphorene nanoribbons[END_REF][START_REF] Yang | Room-temperature magnetism on the zigzag edges of phosphorene nanoribbons[END_REF] makes them particularly intriguing. The highly anisotropic crystal and electronic structure of PNRs opens the possibility for a wide range of spin interactions both in-and out-of-the plane, with their visible band-gap making them well-suited for studying and utilising the interplay between functional optoelectronic properties and magnetism in low-dimensions [START_REF] Ando | Seeking room-temperature ferromagnetic semiconductors[END_REF][START_REF] Gibertini | Magnetic 2D materials and heterostructures[END_REF] .

In this work, we focus on PNRs that are produced readily in batch solutions with typical ribbon widths of ~10 nm and lengths up to 1 μm (see methods) [START_REF] Watts | Production of phosphorene nanoribbons[END_REF][START_REF] Macdonald | Phosphorene Nanoribbon-Augmented Optoelectronics for Enhanced Hole Extraction[END_REF] . The long edges of the PNRs are exclusively zigzag terminated [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] . Figure 1a shows an AFM image of a representative monolayer ribbon (height 0.48 ± 0.14 nm and width 8.12 ± 1.08 nm) along with histograms of the width and height distribution along the ribbon measured using high speed AFM. The green arrow indicates the long direction of the ribbon with zigzag terminated edges. Quantum Monte Carlo and dynamic mean field theory calculations have predicted that such zigzag PNRs host room temperature ferromagnetism on their edges with an energetic difference of just 1 µeV per spin between the edge-to-edge ferromagnetic (FM) and antiferromagnetic (AF) spin aligned ground state [START_REF] Yang | Room-temperature magnetism on the zigzag edges of phosphorene nanoribbons[END_REF][START_REF] Vahedi | Edge magnetic properties of black phosphorene nanoribbons[END_REF][START_REF] Zhu | Magnetism of zigzag edge phosphorene nanoribbons[END_REF] . Furthermore, these edge magnetic states are expected to be robust to external disorder and strain and the widths that we study are in the crossover regime between the decoherence-dominated classical magnetism and fluctuation-dominated quantum magnetism studied in graphene nanoribbons [START_REF] Vahedi | Edge magnetic properties of black phosphorene nanoribbons[END_REF][START_REF] Poljak | Immunity of electronic and transport properties of phosphorene nanoribbons to edge defects[END_REF][START_REF] Golor | Quantum nature of edge magnetism in graphene[END_REF] .

To investigate the magnetic ground state of PNRs, it is important to first establish the degree of intrinsic magnetic anisotropy [START_REF] Sofer | Layered Black Phosphorus: Strongly Anisotropic Magnetic, Electronic, and Electron-Transfer Properties[END_REF][START_REF] Wang | Highly anisotropic and robust excitons in monolayer black phosphorus[END_REF] , due to their highly anisotropic shape and crystal structure. A non-invasive route to accessing the magnetic anisotropy of nanoobjects dispersed in solution is to measure the magnetic field-induced alignment through the linear dichroism (LD) and linear birefringence (LB) [START_REF] Van Rhee | Giant magnetic susceptibility of gold nanorods detected by magnetic alignment[END_REF] . The principle behind optically-detected magnetic alignment is that a magnetically-anisotropic nano-object orients in a magnetic field to minimise its magnetic orientational energy, which can eventually overcome thermal fluctuations. The optical anisotropy (LD or LB) of the underlying nanoobject can then serve as a sensitive probe of the degree of alignment due to the external magnetic field. As shown in Figure 1b, we observe very unexpected behaviour in the magnetic linear birefringence, with the orientational anisotropy saturating at near 10 T at room temperature, much lower than that of other known nanomaterials [START_REF] Pietra | Self-assembled CdSe/CdS nanorod sheets studied in the bulk suspension by magnetic alignment[END_REF][START_REF] Shklyarevskiy | Magnetic Alignment of Self-Assembled Anthracene Organogel Fibers[END_REF][START_REF] Rikken | Manipulation of micro-and nanostructure motion with magnetic fields[END_REF] . On fitting the data to the orientational magnetic energy (see Supplementary Note 1), we retrieve the anisotropy in the volume magnetic susceptibility (∆𝜒) multiplied by the volume of the material, ∆𝜒𝑉 = 5 × 10 !"# 𝑚 $ . We additionally observe a magnetic field induced LD signal at 543 nm, with more light absorbed when the light is polarised along the magnetic field direction (Figure 1b, inset). PNRs are expected to have a higher absorption cross section with light polarised perpendicular to the zigzag direction, as has been reported experimentally for phosphorene and obtained from GW-BSE calculation on PNRs [START_REF] Nourbakhsh | Excitons and optical spectra of phosphorene nanoribbons[END_REF][START_REF] Wang | Highly anisotropic and robust excitons in monolayer black phosphorus[END_REF] . We hence conclude that in a magnetic field, the ribbons orient with their short axis along the field direction.

We then use reported histogram distributions of the length, width and height of 940 PNRs, measured using TEM to calculate the approximate distribution of ribbon volumes in the solution [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] . From this and the calculated ∆𝜒𝑉, we extract the anisotropy in the dimensionless magnetic susceptibility ∆𝜒 ≈ 10 !$ (Figure 1c). This is two orders of magnitude larger than the value expected for a single layer of phosphorene, using extrapolated values from reports of ∆𝜒 for layered phosphorene, i.e., black phosphorous (Figure 1c, green line) . The latter has its predominant magnetic anisotropy between the in-and out-of-plane directions, which are not expected to be the dominant anisotropy directions here [START_REF] Sofer | Layered Black Phosphorus: Strongly Anisotropic Magnetic, Electronic, and Electron-Transfer Properties[END_REF] . We infer that this 100-fold enhancement in the ∆𝜒 must therefore arise from the enhanced number of edges of the ribbons as this is the only intrinsic difference between phosphorene and the PNRs. The origin of the large ∆𝜒 could be from diamagnetic, paramagnetic or ferromagnetic anisotropies, or a combination of these. It is therefore critical to establish the presence of unpaired spins in the PNR system, which would point to a para-or ferromagnetic origin for the large ∆𝜒, and the presence of the predicted edge magnetism.

In order to do this, we performed SQUID magnetometry on a deposited ensemble of drop-cast PNRs in the solid state [START_REF] Buchner | Tutorial: Basic principles, limits of detection, and pitfalls of highly sensitive SQUID magnetometry for nanomagnetism and spintronics[END_REF] . The 7 T-field cooled (FC) and zero field cooled (ZFC) magnetisation measured in a 50 mT field (Figure 1d, inset) shows non-Curie behaviour as we approach temperatures as high as 350 K where the two curves approach each other but do not yet overlap. We performed isothermal magnetic field sweeps at 1.8 and 300 K to study the hysteretic nature of the spin couplings (Figure 1d). We observe an atypical ferromagnetic-like hysteresis at both 300 K and 1.8 K with a coercive field of near 50 mT with the hysteretic behaviour disappearing at small applied fields. This is consistent with the prediction of two competing AF and FM ground states in the system [START_REF] Vahedi | Edge magnetic properties of black phosphorene nanoribbons[END_REF] . Upon applying a small external field of the order of 50 mT, the ferromagnetic interaction energy is lowered by ~ 1 μeV per spin (using 1 Bohr magneton and ∆𝐸 = -𝜇 • 𝐵, for the unpaired electrons in PNRs [START_REF] Zhu | Magnetism of zigzag edge phosphorene nanoribbons[END_REF] ). Hence the degeneracy of the FM and AFM ground state is lifted, and away from zero-field, the FM phase forms the ground state associated with the hysteretic behaviour. This does not change when measured in lower projecting fields of 20 mT (Supplementary Note 2), indicating that the behaviour does not emerge due to an interaction with the measurement/projecting field of the SQUID. In the solution phase (NMP and DMF solvents), we also find that the FC and ZFC magnetisations are split up to the freezing point of the solvent (Supplementary Note 2), above which the liquid dispersions are no longer macroscopic magnetic systems. We repeat the solution measurements on oxygenated PNRs by bubbling air through the solution to chemically terminate the ribbon edge (hydrogenated, phosphate, carboxylate, etc, based terminuses [START_REF] Watts | Production of phosphorene nanoribbons[END_REF][START_REF] Kuntz | Control of Surface and Edge Oxidation on Phosphorene[END_REF][START_REF] Zhu | Stabilizing black phosphorus nanosheets via edge-selective bonding of sacrificial C60 molecules[END_REF][START_REF] Ding | Tuning Phosphorene Nanoribbon Electronic Structure through Edge Oxidization[END_REF] ) and find a purely diamagnetic response (Supplementary Note 2) [START_REF] Peng | Edge effects on the electronic properties of phosphorene nanoribbons[END_REF][START_REF] Hu | Highly Efficient Photocatalytic Water Splitting over Edge-Modified Phosphorene Nanoribbons[END_REF] . Elemental analysis (inductively coupled plasma optical emission spectroscopy (ICP-OES)) of the same PNR solutions estimates less than 20 parts per billion of Fe, Co, Ni, Cr and Mn, ruling these out as the signal origin (Supplementary Note 3). Overall, this suggests that the magnetic response of PNRs arises, as predicted in numerous studies, from the ribbon edge with the persistence of the magnetic behaviour up to room temperature [START_REF] Yang | Room-temperature magnetism on the zigzag edges of phosphorene nanoribbons[END_REF][START_REF] Vahedi | Edge magnetic properties of black phosphorene nanoribbons[END_REF][START_REF] Zhu | Magnetism of zigzag edge phosphorene nanoribbons[END_REF] .

Next, we investigate the local spin correlations along the edge of the PNR by performing continuous wave electron paramagnetic resonance (cwEPR) measurements on a film formed on the inner wall of an EPR sample tube. Consistent with hypothesis of unpaired electrons at the ribbon edge we observe, in Figure 2a,c, a broad ≈12-19 mT wide g≈2 EPR signal at room temperature. To our surprise, a wider magnetic field scan (0-1400 mT) reveals much more intense and even broader FMR like signals (Figure 2a). The EPR and FMR signals both show a strong orientation dependency when rotating around the sample tube long axis as shown in Figure 2b. This behavior arises from the inhomogeneity of the inner wall film in thickness, packing, PNR orientations and size.

The FMR signals seen in Figure 2a,b and the observed orientation behavior is complex due to the inhomogeneity occurring from the inner wall sample. However, in order to shed some light on the internal fields within the sample, we fit the experimental data with two purely Lorentzian and one Voigtian component (see Supplementary Note 4). We find that some PNR domains give rise to net internal fields of ≈ 225-265 mT and are aligned to the external magnetic field direction, thereby shifting the effective resonance position to lower magnetic fields, at around ≈55-95 mT. These PNR domains show very little change in effective resonance positions with changing orientation suggesting they are easily (re-)magnetised. We also see certain PNR domains where the effective resonance position is strongly orientation dependent, ranging from b. Orientational magnetic anisotropy in EPR signal across full range of angles. c. Line cuts of g ≈ 2 signal along with the zero-intensity lines (dashed black) d. Extracted g-factor and linewidth as a function of temperature. We observe a divergence in the linewidth around 70 K characteristic of a phase transition. suggests these domains are hard to (re-)magnetize, resulting in low field signals when the net internal field is (close to) aligned to the external magnetic field direction and high field signals when the net internal field is anti-aligned to the external magnetic field direction. We estimate the net internal field magnitude of these PNR domains to be ≈ 600 mT considering the minimum and maximum resonance positions achieved via sample rotation. PNR domains exhibiting an easy and hard axis could explain the FMR like signals with weak and strong orientating dependency respectively. The broad linewidths of all FMR like signals suggests a wide distribution of internal fields within each domain. This FMR behaviour is strongly reminiscent of reports of magnetism in Ni nanocubes where the different size nanocubes exhibit different internal fields and distributions of internal fields, leading to this rich orientation dependent behaviour [START_REF] Nellutla | Multi-frequency ferromagnetic resonance investigation of nickel nanocubes encapsulated in diamagnetic magnesium oxide matrix[END_REF] (see Supplementary Note 4 for further discussion).

Focusing again on the g≈2 EPR signal, we see that the g-factor and linewidth vary as a function of sample tube rotation, with the g-factor ranging from 2.03 to 2.21 (see Supplementary Note 4). The g-factor for black phosphorous has been calculated to be 2.0 in-plane and 2.14 out-ofplane in the monolayer (1L) limit, with a strong layer number dependence, reaching 2.9 in the bulk out-of-plane limit [START_REF] Zhou | Effective g factor in black phosphorus thin films[END_REF] . The calculated out-of-plane g-factor for a 2 layer (2L) black phosphorous is ≈ 2.21. Hence the observed EPR signal, results predominantly from 1L and 2L PNRs, consistent with Watts et al. [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] , and the spin polarisation is either in-plane or out-of-plane depending on the tube orientation.

The EPR signal also shows a strong temperature dependency, where the linewidth and g-factor behaviour are non-monotonic, with the divergence of the linewidth and the disappearance of the signal around 55 K, which reappears at lower temperatures. This behaviour which also appears, albeit weakly, in the SQUID magnetisation (Figure 1d, inset), is characteristic of a phase transition, however further work is required to describe its nature. We note that the Raman spectrum shows no significant evolution with temperature (300 K to 10 K) that would be indicative of a structural phase transition (Supplementary Note 5). We suggest that the smooth temperature dependence of the g-factor from 300 to 75 K could signal a slow spincanting in the out-of-plane direction due to AF correlations which lowers the g-factor from the fully out-of-plane polarised value of 2.14 towards the in-plane value of 2.00 [START_REF] Mills | Surface anisotropy and surface spin canting in the semi-infinite ferromagnet[END_REF] . There appears to be a broad photobleaching spectra signature at 700 nm that overlaps with the ESA bands. b. Polarisation and magnetic field resolved decay of ESA 1 kinetics. Upon application of a 800 mT field, the ribbons align in the magnetic field, splitting the degeneracy of the HH and VV polarised experiments. Further, we find a significant reduction in the signal magnitude when photoexciting V polarised (along the PNR long axis). c. Magnetic field induced anisotropy obtained by subtracting the spectral response of the unaligned ribbons. This unveils a spectral signature centred at 850 nm (corresponding to ESA 1), the magnitude of which can be strongly modulated with a sub-1T field (right panel). d. Decay of polarisation anisotropy (r) as a function of magnetic field. In an unpolarised solution under no magnetic field, we find that the anisotropy is initially -0.2 and decays to close to 0 within 1 ps (green and blue), indicating the photoexcited state is perpendicular to the probed state after 1 ps (the differences in the V and H polarised experiments at 0 mT demonstrate the setup's resolution). Upon application of the 800 mT field however, photoexciting V polarised (PNR long axis, orange curve) we find that the anisotropy retains its' sign (indicating the photoexcited state is the same orientation as the probed one) whereas when photoexciting H polarised (PNR short axis, red curve), the anisotropy flips sign, indicating a dipole flip, as in the unpolarised case, with the average converging to the unpolarised system, as expected.

Having established the signatures of key spin correlations and room temperature magnetism in PNRs, we turn to their optical properties. The temperature dependent absorption spectrum of PNRs is broad and featureless (see Supplementary Note 6), we hence interrogate the optical properties by studying the photophysical dynamics [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] . Upon photoexcitation with a 460 nm, 250-fs pulse, we see the emergence of two long lived (nanosecond) excited state absorptive (ESA) bands at 800 nm (ESA 1) and 500 nm (ESA 2) along with an overlapping, short-lived (picosecond) broad photobleaching band at 700 nm (Figure 3a).

From the LD in a magnetic field (Figure 1b, inset), we know that the ribbons align with their short axis along the magnetic field direction in solution. In Figure 3b, we show the results of performing polarization resolved pump-probe experiment on the PNRs in solution, with an 800 mT applied external magnetic field in the Faraday geometry. We find that the typical copolarised horizontal (HH) and vertical pump and probe (VV) signal's degeneracy is lifted and a lab axis is introduced into the experiment upon application of an 800 mT field. The magnetic field is polarised along the H direction in the lab axis which, as expected, results in the higher transient absorption signal due to the increased net linear dichroism along the HH direction (Figure 1b). Further, as expected, we see the alignment effect scales with the applied field strength (Supplementary Note 7).

We leverage this magnetic alignment effect to isolate the spectral features of the magnetically aligned ribbons (as opposed to optical signatures of a non-magnetic sub-population) by subtracting the HH and VV polarised transient absorption maps at different field strengths, as any unaligned ribbons would have identical HH and VV polarised transient absorption maps yielding no difference signal. As seen in Figure 3c, we find that the magnitude and kinetics of the photoinduced absorption band (ESA 1) is most strongly perturbed by the external field. We can therefore confidently assign ESA 1 to that of the magnetically active ribbons.

Having established a spectral signature of the magnetically active ribbons, we use this feature to track the evolution of the photoexcitation in time in the magnetically active ribbons through the transient absorption anisotropy in the ESA 1 band (800-900 nm) defined as, [START_REF] Min | Transient absorption anisotropy study of ultrafast energy transfer in porphyrin monomer, its direct meso-meso coupled dimer and trimer[END_REF] . We find that at 200 fs the anisotropy is -0.2, signifying that we are likely probing the dipole polarised orthogonally to the main photoexcited transition. Over 1 ps, the anisotropy relaxes to a value of about 0.05, signifying that there is energetic transfer between the dipoles as the timescale is faster than the ns-orientational relaxation dynamics typically associated with molecular anisotropy. As expected, we demonstrate the H and V pumped anisotropies are degenerate, demonstrating the instrument resolution (Figure 3d, blue and green).

𝑟 = ∆" " ∥ ! ∆" " $ ∆" " ∥ &" ∆" " $
Upon application of the 800 mT field due to the macroscopic alignment of the ribbons, we are able to increase our dipole selectively and pump the short-axis polarised (H) and long-axis polarised (V) dipole. It is worth noting that as the ribbons are now aligned, the sign of the anisotropy is the key feature of interest, as the value of the anisotropy will not be consistent with the isotropic molecular distribution assumption standard to such measurements 57 . We find that when pumping the long-axis polarised (V) dipole (Figure 3d, orange), the anisotropy does not change sign, indicating that the energy remains in this state, whereas when we pump the short-axis polarised (H) dipole, the sign flips (red), just like in the case of the unpolarised ribbons (blue). Further, the fact that we can construct the unaligned PNR anisotropy as an average of pumping the 800 mT aligned long-axis polarised (V) and short-axis polarised (H) dipoles (Figure 3d, purple) confirms that these are the main two dipoles involved in the optical response of the system. Taken together with the larger magnitude of the transient absorption signal (Figure 3b) of the H pumped system, this suggests that the photoexcitation is initially predominantly polarised along the short axis and transfers to an orthogonal dipole polarised along the PNRs long axis on a sub-1ps timescale. shows that, as previously reported [START_REF] Ribeiro | Edge phonons in black phosphorus[END_REF] , the B3g 1 mode is localised on the edge. b. Wavelength resolved, impulsive, time-domain Raman spectrum of the photoexcited PNRs reveals that the B3g 1 mode is most strongly coupled (white curve) to ESA 1 (red curve) which was also most strongly modulated by the external magnetic field.

GW-BSE calculations of the optical spectrum of zigzag PNRs report that the lowest excited state is a dark exciton which is polarised along the ribbons long axis [START_REF] Nourbakhsh | Excitons and optical spectra of phosphorene nanoribbons[END_REF] . We propose that the primary photoexcitation mechanism is therefore that the initial excitation dipole polarised along the ribbon short axis decays to this long axis polarised dark exciton on sub-1ps timescales. This is consistent with the fast disappearance of the photo bleach signal corresponding to the initial state in the transient absorption (Figure 3a) with the formation of long lived ESAs likely from the dark exciton.

While calculations have reported that the first bright exciton in zigzag PNRs has significant density along the ribbon edges, it is unclear whether this lower energy dark exciton that the photoexcitation decays to is also localised along the ribbon edges [START_REF] Nourbakhsh | Excitons and optical spectra of phosphorene nanoribbons[END_REF] . This is critical for any potential excited state magneto-optical coupling in PNRs, as most of the photoexcitation lifetime appears to be spent in this state. We exploit the presence of special localised phonon modes on the edges of black phosphorous to study this [START_REF] Ribeiro | Edge phonons in black phosphorus[END_REF] .

The spontaneous Raman spectrum of the PNRs (Figure 4a, blue curve) is dominated by the A 1 g, B 2 g and A 2 g modes, all of which are also observed in bulk black phosphorous (BP) [START_REF] Guo | From Black Phosphorus to Phosphorene: Basic Solvent Exfoliation, Evolution of Raman Scattering, and Applications to Ultrafast Photonics[END_REF] . We demonstrate the previously-reported, bulk-symmetry forbidden B 1 g and B3g 1 edge phonon modes along the zigzag edge of BP using Raman imaging (Figure 4a, inset) [START_REF] Ribeiro | Edge phonons in black phosphorus[END_REF] . In order to study the mode coupling to the excited state energy surface, we perform impulsive vibrational spectroscopy (IVS) using broadband compressed 13-fs pulses centred at 550 nm. We find that when resonantly photoexciting the PNRs, the photoexcited state is strongly coupled to the B3g 1 edge phonon mode (Figure 4a, red curve). Control experiments with off-resonant IVS on the PNRs and DMF solvent confirm that the B3g 1 edge phonon mode is indeed dominantly coupled to the excited state potential energy surface (see Supplementary Note 8). Furthermore, as seen in Figure 4b, we find that this edge phonon mode is most strongly coupled to ESA 1, the transition characteristic of the magnetically active ribbons. The mode selectivity of the excited state coupling to edge mode as compared to the bulk modes suggests that the electron density of the photoexcitation is along the edges of the PNRs. In the solution phase, we observe that the excited state IVS of the B3g 1 out-of-plane edge mode is present when studying the unaligned (0 mT) ribbon solution only in the perpendicular pump-probe configuration, i.e. HV or VH (Supplementary Note 8) and not the parallel configurations. In other words -one needs to probe a dipole perpendicular to the main photoexcitation (that we ascribe to this dark excitonic transition polarised along the PNRs long axis) in order to observe the excited state moving along the B 3 g 1 vibrational coordinates at the edge of the PNRs.

In summary, we have evidenced via SQUID magnetometry and cwEPR, macroscopic magnetism at room temperature in films of PNRs arising from the ribbon edges. In solution a giant magnetic anisotropy means PNRs can be readily aligned with their short axis along the field direction. On photoexcitation, energy rapidly migrates from a bright-to dark-exciton state that is localised to the PNR edge and strongly coupled to a symmetry-forbidden edge phonon mode. Our findings open up extensive opportunities and challenges for further understanding the magnetic properties of PNRs (spin coherence times, influence of ribbon width/strain, role of thickness, etc) and for the building of proposed nanoribbon devices [START_REF] Macdonald | Phosphorene Nanoribbon-Augmented Optoelectronics for Enhanced Hole Extraction[END_REF][START_REF] Zhao | Computational Study of Tunneling Transistor Based on Graphene Nanoribbon[END_REF][START_REF] Baringhaus | Exceptional ballistic transport in epitaxial graphene nanoribbons[END_REF][START_REF] Wang | Graphene nanoribbons for quantum electronics[END_REF][START_REF] Son | High-Performance Broadband Faraday Rotation Spectroscopy of 2D Materials and Thin Magnetic Films[END_REF] as well as more disparate applications (non-linear optical elements 64-66 , nanorobotics [START_REF] Rikken | Manipulation of micro-and nanostructure motion with magnetic fields[END_REF][START_REF] Zeng | Autonomous Catalytic Nanomotors Based on 2D Magnetic Nanoplates[END_REF][START_REF] Nguyen | Magnetic control of graphitic microparticles in aqueous solutions[END_REF] , magnetic drug delivery [START_REF] Shaw | Antipathogenic properties and applications of low-dimensional materials[END_REF] ). Important questions are raised about the exact microscopic arrangement of spins at the ribbon edge, the in/out-of-plane couplings and couplings between layers. More generally, it remains to further directly link the semiconducting and magnetic properties of PNRs, but the optically active nature of PNR edges provides a promising way for doing so via excitation of either the electronic state or edge phonons.

Methods

PNR Synthesis

PNRs were synthesised as reported in Watts et al. [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] . In short, in an argon glovebox (<0.1 ppm Oe, <0.1 ppm H2O) black phosphorus (2D semiconductors Ltd and Smart Elements) was ground via pestle and mortar into ca. 1 mm flakes and 124 mg was transferred to a glass tube fitted with a metal Swagelok valve, alongside 3.5 mg freshly cut lithium metal (Sigma Aldrich, SA, 99% rod). The tube was evacuated (ca. 10 -7 mbar) and cooled to -50 o C and ammonia gas (SA, 99.95%, precleaned by condensation over excess lithium metal) was condensed to submerge the lithium and phosphorus. The solution immediately turned dark blue from the formation of lithium electride, and slowly turned orange over 16 h. The ammonia was then evaporated and the orange LiP8 salt was dried under vacuum (ca. 10 -7 mbar) at room temperature. The LiP8 (10 mg) was then placed in 10 mL of either NMP (SA, 99.5% anhydrous) or DMF (SA, 99.9% anhydrous) which had been dried with 4 Å molecular sieves for 1 week. The mixtures were bath sonicated for 30 min, before centrifuging (100g, 10 min) and decanting in a glovebox to give solutions of PNRs.

Atomic Force Microscopy

PNR samples were prepared via drop-casting from solution onto freshly cleaved graphite (HOPG) substrates as discussed previously [START_REF] Watts | Production of phosphorene nanoribbons[END_REF] . High-resolution topography maps were then collected using contact-mode HS-AFM (Bristol Nano Dynamics Ltd, UK) with silicon nitride microcantilevers (MSNL-10, Bruker) with nominal tip radii of 2 nm. To ensure that the images obtained featured tip-limited lateral resolution scan sizes of 0.8 µm × 0.5 µm were used, corresponding to pixel sizes of ~1 nm × 1 nm. HS-AFM images containing ribbons selected for analysis were input into custom LabVIEW image analysis software (National Instruments) that isolated the ribbon from the background in an automated fashion before generating histograms of ribbon height and width. In brief, the algorithm consisted of: a local thresholding step to identify the ribbon; a mask step to separate pixels corresponding to the background from the ribbon; an erosion step to remove the edge pixels where the AFM tip traverses the slope between the background and ribbon; calculation of height histograms for both the background and ribbon pixels and finally; Gaussian fits to the histograms recovered the mean height of each and enabled the calculation of the height of the ribbon in respect to the background.

The width of the ribbon was analysed in parallel using the same algorithm. After the masking step the central axis of the isolated ribbon was determined by fitting cubic splines to the two edges of the ribbon and calculating the mean of these two edges. Then, the width of the ribbon was determined by calculating the distance to the central axis from each edge. Measurements of the width of the ribbon along its entire length were then input into a histogram and fit with a Gaussian distribution to identify the mean width of the PNR.

Magnetic Linear Dichroism and Birefringence

The degree of magnetic alignment is measured optically, through the magnetic field-induced linear dichroism (LD) and birefringence (LB), using red and green Helium Neon lasers (wavelengths 632.5 and 543 nm, respectively). The solution samples were contained in an optical cuvette (thickness 5 mm or 10 mm) positioned within a temperature-controlled environment at 20.0 ± 0.1 °C in a 33 T Florida-Bitter electromagnet or a Varian V-3900 2 T magnet. The LD and LB signals were measured using standard polarization modulation techniques using a photo-elastic modulator [START_REF] Kemp | Polarized Light and its Interaction with Modulating Devices-a Methodology Review-Polarized Light and Its Interaction With Modulating Devices[END_REF] .

SQUID Magnetometry

Magnetisation measurements were obtained in a Quantum Design Magnetic Properties Measurement System (MPMS 3) using a Superconducting Quantum Interference Device (SQUID) magnetometer. Measurements were performed down to 1.8 K and up to 350 K, in various applied magnetic field strengths up to 7 T. Data for two isothermal magnetic field scans of PNRs out of solution are shown in Figure 1d, and two magnetisation-temperature curves of PNRs are displayed in the inset of Figure 1d. Each magnetisation-temperature curve was measured starting at 1.8 K and increasing the temperature in a constant probing field of 50 mT. Prior to beginning any new magnetic field or temperature scan the system was taken to 300 K and the magnet was reset to remove any stray flux from the SQUID. Samples were secured on an MPMS quartz sample holder using GE varnish, and care was taken to ensure that they were not touched by any magnetic material throughout the mounting and loading procedure. At each measurement, several dc magnetisation measurements were averaged, providing a reliable measure of the bulk magnetisation of the PNR samples.

While measuring the PNRs in the DMF and NMP solutions (Supplementary Note 2) 10 μL were pipetted into a plastic capsule. The capsule was fixed to a quartz rod using GE varnish, and the presence of the PNR increased the magnetisation by two orders of magnitude when compared to the GE and capsule by itself. The dc moment was averaged over multiple scans, and between each experimental run the MPMS was brought to 300 K and the magnet was reset to remove any remaining flux in the SQUID. Oxygenation of the PNR was performed by introducing by bubbling an O2 gas stream to the PNR in the plastic capsule while ensuring it does not come into contact with any metallic substances.

Electron Paramagnetic Resonance

For the cwEPR experiments, 100 µL of PNR in DMF solution was filled into a 4 mm outer diameter (3 mm inner diameter) quartz EPR tube, inside a nitrogen glovebox. The sample tube was attached to a custom adaptor and transferred to a pumping station outside the glovebox. The custom adaptor keeps the sample in the inert glovebox environment. The solution was evaporated under vacuum (using a pumping station), resulting in a film on the inner walls of the EPR tube. The procedure was repeated four times, to achieve a thicker film, in total 400 µL was evaporated. The inner wall sample was then left to pump to a pressure of 6×10 -4 mbar and flame-sealed.

Narrow Magnetic Field Scan -Rotation and Temperature Series (Setup 1):

The cwEPR spectra were recorded at X-band (≈9.4 GHz) using a laboratory-built EPR spectrometer (all narrow scan cwEPR spectra were rescaled to 9.4 GHz). The setup consists of a Bruker ER 041 MR microwave bridge together with an ER 048 R microwave controller and an AEG electromagnet together with a Bruker BH 15 Hall effect field controller. The magnetic field was additionally monitored with a Bruker ER 035 M NMR Gaussmeter. The resonator used was a Bruker ER 4122-SHQE resonator. The static magnetic field was modulated at 100 kHz and lock-in detection was carried out using a Stanford Research SR810 lock-in amplifier in combination with a Wangine WPA-120 audio amplifier. An ESR 900 helium flow cryostat together with a ITC503 temperature controller (Oxford Instruments, UK) was used for low temperature measurements. The spectra were acquired at a frequency of ≈ 9.4 GHz with a microwave power of 7.96 mW and 1 mT modulation amplitude. The magnetic field was calibrated using a standard N@C60 sample with a known g-factor.

The narrow magnetic field scans presented in the main text (Figure 2c) have had the FMR like signal (slope) removed from them using an appropriate polynomial fit. See Supplementary Note 4 for further details.

Wide Magnetic Field Scan -Rotation Study (Setup 2):

The cwEPR spectra were recorded at X-band (≈9.55 GHz) using a laboratory-built EPR spectrometer (All wide scan cwEPR spectra were rescaled to 9.55 GHz). The setup consists of a Bruker ER 046 XK -T microwave bridge together with an ER 048 R microwave controller and a Varian electromagnet together with a Bruker ER 032M Hall effect field controller. The resonator used was a Bruker MD5 dielectric ring resonator. The static magnetic field was modulated at 99 kHz and lock-in detection was carried out using a Stanford Research SR830 lock-in amplifier in combination with a Wavetek 50 MHz function generator model 80. The wide magnetic field measurement was carried out at 296 K. The spectra were acquired at a frequency of ≈ 9.55 GHz with a microwave power of 7.96 μW and 0.4 mT modulation amplitude.

Transient Absorption Spectroscopy

The transient absorption measurements from 550 -930nm were performed using a homebuilt setup around a Yb:KGW amplifier laser (1030 nm, 38 kHz, 15 W, Pharos, LightConversion). The probe pulse was a chirped seeded white light continuum created using a 4 mm YAG crystal that spanned from 500 nm to 950 nm. For the source of the pump pulse (∼200 fs) a commercial optical parametric amplifier (OPA) OPHEUS ONE (Light Conversion) was used.

The transient absorption measurements from 400 -550 nm were performed using a homebuilt setup around a Ti-Sapphire (800 nm, 1 kHz, Spectra-Physics, Solstice Ace). The probe pulse was a chirped seeded white light continuum created using a 4 mm CaF2 crystal that spanned from 400 nm to 600 nm. For the source of the pump pulse (∼100 fs) the fundamental of the laser was doubled in a β-barium borate crystal.

The sub-1T magnetic field was generated using an electromagnet from GMW Model 3470 with 1 cm distance between cylindrical poles and the field strength calibrated with a Gaussmeter.

Impulsive Raman Spectroscopy

A femtosecond time domain Raman spectroscopy measurements were performed using a homebuilt setup around a Yb:KGW amplifier laser (1030 nm, 38 kHz, 15 W, Pharos, LightConversion). The probe pulse was a chirped seeded white light continuum created using a 4 mm YAG crystal that spanned from 500 nm to 950 nm. The pump pulse for the resonant experiment was created using a non-collinear optical parametric amplifier (NOPA) where the 1030 nm seeded a white light continuum stage in sapphire which was subsequently amplified with the third harmonic of the 1030 nm laser in a β-barium borate crystal to create a broad pulse centred at 550 nm. The pump pulse for the off-resonant experiment was created using a non-collinear optical parametric amplifier (NOPA) where the 1030 nm seeded a white light continuum stage in a YAG crystal which was subsequently amplified with the third harmonic of the 1030 nm laser in a β-barium borate crystal to create a broad pulse centred at 750 nm. Both pulses were compressed using a chirped mirror and wedge prism (Layerterc) combination to a temporal duration of under 15 fs. Compression was determined by second-harmonic generation frequency-resolved optical gating (SHG-FROG; upper limit) and further confirmed by reference measurements on acetonitrile where the 2,200 cm -1 mode could be resolved. The probe white light was delayed using a computer-controlled piezoelectric translation stage (Physik Instrumente), and a sequence of probe pulses with and without pump was generated using a chopper wheel (Thorlabs) on the pump beam. The average fluence of the pump 10 μJ /cm 2 .

Temperature dependent Raman spectroscopy and Raman Imaging

Raman spectra were measured as a function of temperature from 4 K to 300 K. Raman measurements were conducted by backscattering (T64000, Horiba) a CW diode line (532 nm, 1 mW). Spectra were collected at > 200 cm -1 , where the CCD detector (Horiba Synapse Open-Electrode) has a monotonically increasing quantum efficiency of 0.43 -0.50. Acquisitions employed a 100× optical objective and used minimal laser intensity to avoid sample degradation.

For Raman imaging a standard layout of an epi-detected Raman microscope was used. A pump laser beam (wavelength = 532 nm, Coherent Mira) was spectrally cleaned up by a bandpass filter (FLH05532-4, Thorlabs), and its beam width was expanded to 7.2 mm before entering a home-built inverted microscope. Additional waveplates (half-waveplate and quarter-waveplate for 532 nm, Foctek Photonics) pre-compensated the ellipticity introduced by the dichroic filter (F38-532_T1, AHF) and also generated circularly polarised light. We used high numerical aperture (NA) oil-immersion objectives (Nikon 60X/1.4 N.A. oil) to ensure high-resolution imaging and increase collection efficiency. Pump power before the objective was 30 mW, a power level that ensured no degradation of samples. The samples were scanned with galvanometric mirrors (Thorlabs). The Raman inelastic backscattered light was collected by the same objective and focused with the microscope tube lens either onto the slit of a spectrometer (Andor, Shamrock 303i, grating 300 l/mm; the slit also acts an effective pinhole for confocal detection). The spectrometer is equipped with a high-sensitivity charge-coupled camera (Andor, iXon 897). All images presented were taken with integration times/pixel in the 500 ms range. Recording of data was performed by a custom Matlab program. For Raman imaging, black phosphorous flakes (300-600 nm thickness; 2D semiconductors Ltd) were mechanically exfoliated inside a nitrogen glovebox and transferred onto a Si substrate. A 0.15 mm thick coverslip was placed over the flakes and sealed with epoxy glue to act as an encapsulant. The polarisation of the pump and Raman light was not strongly controlled, but simply adjusted to maximise the respective signals.

Figure 1 :

 1 Figure 1: Shape, magnetic and optical anisotropy of PNRs: a. (Top) Crystal structure of the synthesised PNRs with the zigzag edges along their long axis (green arrow). (Bottom) Representative AFM image and histogram of height and width along a single ribbon. An average width of 8 nm and average monolayer height of 0.6 nm is found along the ribbon. b. Alignment of PNRs in solution measured using magnetic field induced linear birefringence up to 30 T shows a saturation at 10 T. Combined with the magnetic linear dichroism (inset), these result suggest that PNRs align with their short axis along the field direction. c. Calculated histogram of the intrinsic anisotropy in the magnetic susceptibility shows a 100-fold edge enhancement from the previously reported monolayer phosphorene limit. d. Isothermal magnetisation of PNR at 1.8 K (blue) and 300 K (red) measured using a SQUID magnetometry. Ferromagnetic-like correlations up to room temperature are demonstrated by the hysteresis. The inset displays the field-cooled (FC) at 7 T and zero-field cooled (ZFC) magnetisation as a function of temperature for the same PNR samples.

Figure 2 :

 2 Figure 2: EPR and FMR-like local magnetic correlations: a. cwEPR signal of PNRs in the solid state measured over a wide magnetic field range shows both the unpaired electron peak at 320 mT sitting atop of a much larger ferromagnetic-like resonance at 82 mT, suggesting an internal field of 240 mT (green curve; 15 degrees -low field signal: internal field aligned to magnetic field). The sample is highly inhomogeneous hence rotation results in a shift in the resonance positions (red curve; 105 degrees -high field signal: internal field anti-aligned to magnetic field). b. Orientational magnetic anisotropy in EPR signal across full range of angles. c. Line cuts of g ≈ 2 signal along with the zero-intensity lines (dashed black) d. Extracted g-factor and linewidth as a function of temperature. We observe a divergence in the linewidth around 70 K characteristic of a phase transition.

Figure 3 :

 3 Figure 3: Photoexcited dipole relaxation dynamics. a. Broadband transient absorption spectrum of PNRs. Upon photoexcitation the PNRs display two broad excited state absorption (ESA) bands that decay on the order of ns.There appears to be a broad photobleaching spectra signature at 700 nm that overlaps with the ESA bands. b. Polarisation and magnetic field resolved decay of ESA 1 kinetics. Upon application of a 800 mT field, the ribbons align in the magnetic field, splitting the degeneracy of the HH and VV polarised experiments. Further, we find a significant reduction in the signal magnitude when photoexciting V polarised (along the PNR long axis). c. Magnetic field induced anisotropy obtained by subtracting the spectral response of the unaligned ribbons. This unveils a spectral signature centred at 850 nm (corresponding to ESA 1), the magnitude of which can be strongly modulated with a sub-1T field (right panel). d. Decay of polarisation anisotropy (r) as a function of magnetic field. In an unpolarised solution under no magnetic field, we find that the anisotropy is initially -0.2 and decays to close to 0 within 1 ps (green and blue), indicating the photoexcited state is perpendicular to the probed state after 1 ps (the differences in the V and H polarised experiments at 0 mT demonstrate the setup's resolution). Upon application of the 800 mT field however, photoexciting V polarised (PNR long axis, orange curve) we find that the anisotropy retains its' sign (indicating the photoexcited state is the same orientation as the probed one) whereas when photoexciting H polarised (PNR short axis, red curve), the anisotropy flips sign, indicating a dipole flip, as in the unpolarised case, with the average converging to the unpolarised system, as expected.

  absorption signal with the probe polarisation parallel and perpendicular to the pump polarisation

Figure 4 :

 4 Figure 4: Coupling of the photoexcitation to the magnetic edges a. Excited state (impulsive, time domain) Raman spectrum of the PNRs (red) shows a strong coupling to the B3g 1 mode compared to the ground state spectrum (blue), indicating that the photoexcitation is delocalised along the edges of the ribbon. Inset: (Left) Zoom into low frequency region of ground state Raman spectrum. (Right) Raman imaging of a black phosphorous flakeshows that, as previously reported[START_REF] Ribeiro | Edge phonons in black phosphorus[END_REF] , the B3g 1 mode is localised on the edge. b. Wavelength resolved, impulsive, time-domain Raman spectrum of the photoexcited PNRs reveals that the B3g 1 mode is most strongly coupled (white curve) to ESA 1 (red curve) which was also most strongly modulated by the external magnetic field.

Funding: A.A. acknowledges funding from the Gates Cambridge Trust as well as support from the Winton Programme for the Physics of Sustainability. S.F. acknowledges funding from the Studienstiftung des deutschen Volkes and the Engineering and Physical Sciences Research Council (EPSRC UK) via an EPSRC Doctoral Prize Fellowship. R.P. acknowledges Clare College, Cambridge for funding via a Junior Research Fellowship. We acknowledge financial support from the EPSRC via grants EP/M006360/1 and EP/W017091/1 and the Winton Program for the Physics of Sustainability. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 758826). This work was supported by HFML-RU/NWO-I, member(s) of the European Magnetic Field Laboratory (EMFL) and by EPSRC (UK) via its membership to the EMFL (grant no. EP/N01085X/1).

Data and materials availability:

The data that support the plots within this paper and other findings of this study are available at the University of Cambridge Repository (https://doi.org/XXXXX).

Competing interests:

The authors declare that they have no competing interests.