text:
# Hard X-Ray Lightcurves of High Mass X-Ray Binaries

## 1 Introduction

The Burst And Transient Source Experiment on the Compton Gamma Ray Observatory \[Fishman et al 1989\] provided near continuous monitoring of the whole sky at hard X-ray energies. BATSE only sees the high energy tail of HMXBs, typically detected in the lowest DISCLA channel and CONT channels 1-4 covering 20-70keV. In studying X-ray pulsars two modes of operation are possible; pulsed flux and Earth occultation. If pulsations are detected corresponding to the spin period of a known pulsar then the pulsed flux and pulse frequency can be measured with Fourier analysis or epoch folding (Bildsten, L. et al. 1997) and used to monitor the source with great sensitivity and fine time resolution. \[Bildsten et al 1997\] For the long term study of behavior related to orbital motion in binary systems, daily averages of the total flux are sufficient and are in addition sensitive to non-pulsed X-ray output from the source. By exploiting the Earth’s limb as an occulting mask, individual X-ray sources can be isolated with an accuracy of 1$`\mathrm{°}`$. Termed the Earth Occultation Technique, two measurements of the total flux from a source are made during each 90 minute satellite orbit as the source moves into and out of eclipse due to the satellite’s motion, full details are given by Zhang et al \[Zhang et al 1994\] & Harmon et al \[Harmon et al 2002\].

In this paper we present folded X-ray lightcurves of Type I (normal) outbursts from 8 HMXBs in which periodic fluctuations corresponding to previously published orbital periods were significantly detected. Such outbursts are thought to occur at each periastron passage of the neutron star. These outbursts last for 8-10 days and the X-ray luminosity increases by a factor of $``$10. In contrast, Type II are giant outbursts which do not correlate with the orbital phase but tend to last for several orbital cycles. The X-ray luminosity during the outburst may reach close to the Eddington luminosity.

In addition we studied the 30.5d modulation from LMC X-4 which was also a significant detection and investigated the X-ray behaviour of A0538-668 in the light of recent optical work \[Alcock et al 1999\]. The BATSE datasets for these systems cover many orbital cycles and hence present the first ever opportunity to study their long term behavior. Previously only widely spaced, scheduled or chance obsevations were available for this purpose. Usually made with different instruments having little or no overlap in terms of spectral response, extrapolating long term trends from previous data was of limited value.

The publicly available pulsed flux monitor reports produced by MSFC demonstrates that HMXBs typically exhibit active epochs lasting several months interspersed with quiescent periods. A question we addressed was that of possible changes in behaviour between these epochs. If a source is bright enough or its active epochs span enough orbital cycles then folded X-ray lightcurves from different epochs can be directly compared and variations in behavior detected.

Variability in the X-ray recurrence period can be looked for in the power spectrum of the flux history. Width of the recurrence frequency peak will be increased if that frequency varies, in other words if the exact time at which the outbursts turn on varies from orbit to orbit then this will show up in the power spectrum. If such jitter can be identified this could provide clues about the underlying physical behaviour of such systems. For example theory indicates that accretion rates are effected by a neutron star’s speed relative to its surroundings. The Alfven speed in stellar wind plasma is a function of both magnetic field strength and plasma density. Thus a high outflow velocity of a wind or circumstellar disk may inhibit accretion during periastron passage and push the X-ray outburst to a later phase. Additionally the orientation of the magnetic field frozen into a plasma strongly affects its coupling to a magnetosphere moving through it. These effects, although extremely difficult to measure directly, may leave a signature in possible chaotic fluctuations or jitter in the exact phase at which outbursts turn on. The jitter is expected to be chaotic because there is no direct connection between the plasma properties and the binary motion, yet the stability of the accretion disk depends on a dynamic balance between accretion rate onto the disk and accretion rate onto the neutron star.

Secondary outbursts occuring close to apastron of the neutron star are also a key observation in this study. Their existence and exact phase at which they occur will yield information about circumstellar disk outflow and stellar wind parameters, principally outflow velocity and density information. Higher outflow velocity is expected to delay periastron outburst turn-on and is cited as a possible cause of apastron outburst in the Bondi-Hoyle theory of accretion. Our folded lightcurves show that EXO2030+375, GX301-2 and 4U1145-619 may exhibit X-ray emision at apastron passage of the neutron star. Whether this emission occurs always or only sporadically, is impossible to determine from individual orbits, so long-term coverage is a necessity.

In these preliminary interpretations of our lightcurves we are aware that BATSE only sees the high energy tails of the emitted X-ray spectra of these sources. Lightcurve shape is thus affected by the true spectral profile of the emitted X-rays and the spectral response of the detector. The most significant effect is expected to occur at the peak of periastron outbursts, when there is the highest column density of obscuring matter and the most intense flux of hard photons, leading to X-ray reprocessing and re-radiation at lower energies.

Parameters for the 8 systems studied are shown in Table 1. Of these sources, in all cases but LMC X-4 the modulation studied is directly linked to orbital motion. In LMC X-4 the binary period of 1.4d is not seen because it is too faint in hard X-rays. The 1.4 day orbital period of LMC X-4 is above the Nyquist limit for this data (about 2.4d) so aliasing will shift this frequency in the power spectrum, but no significant signal is detected. However the well known periodicity near 30.5 days was detected and a lightcurve constructed. An improved measurement was made of X-ray modulation due to the precession of a warped accretion disk.

Note that throughout this paper dates are presented as Truncated Julian Dates (TJD), where TJD = JD - 2440000.5.

## 2 Data processing and data artefacts

The BATSE dataset used in this work consists of daily averages of the flux levels derived from the occultation technique at MSFC \[Harmon et al 2002\].

In the occultation technique, source fluxes are in principle measured twice every spacecraft orbit, as rising and setting steps. However, in practice, not all occultation steps are available for measurement. The most common reason steps are unavailable are (1) passage of the spacecraft through the lower Van Allen radiation belt at the South Atlantic Anomaly, when the detector high voltage was turned off; (2) periods of time when CGRO was out of line-of-sight contact with the NASA Tracking and Relay Satellites (TDRS); and (3) data which were flagged by the BATSE operations team (i.e. during gamma ray bursts, solar flares, etc.) and were not available for analysis. The oblateness of the Earth causes the CGRO orbit to precess with a period of about 53 days. During the precession cycle, high declination sources ($`|\delta |41^{}`$) cease occulting near the orbital poles.

When occultation steps are measured, the background is estimated using 110 seconds of data before and after the occultation step, called the fitting window. Occultation steps from sources other than the measured source within this window are considered interfering sources. If occultation steps of an interfering source occur within 10 seconds of the measured source, that step is discarded. The number of interfering sources in the fit affects the systematic error in three ways (1) where sharing of flux occurs between source terms in the fitting process and (2) where sources are either included or not included in the fit because of an incorrect assumption about the sources’ relative intensity to the measured source or (3) when an unknown source is present in the fitting window. Other sources of systematic error include red noise and coherent pulses from bright X-ray sources present, but not necessarily occulting, within the fitting window.

To estimate the average systematic error, a grid covering the galactic plane with points every 3 degrees along the galactic plane with two additional sets of grid points at +6 degrees and -6 degrees galactic latitude. Grid points within 2 degrees of bright occultation sources where excluded, yielding a total sample of 156 points. The average flux from each grid point was averaged over a 9 year period and the standard deviation of 1-day flux averages was computed as a function of galactic longitude. The values for the standard deviations are larger than expected for a Gaussian distribution about zero flux due to various systematic effects. We find broadening factors of about 30%-60% in excess of normal statistics, strongest near the galactic center and near bright variable sources such as Vela X-1 and Cygnus X-1. Both positive and negative trends are present in the average fluxes. For example, near EXO2030+375 ($`l_{II}`$ = 77.2 degrees, $`b_{II}`$ = -1.3 degrees), the average flux at $`b_{II}`$ = 0 degrees is 0.0077, 0.0009, -0.001 photons cm$`{}_{}{}^{2}s_{}^{1}`$ and standard deviation 1.85, 1.42, 1.26 (1 sigma = 0.01 photons cm$`{}_{}{}^{2}s_{}^{1}`$) for longitudes $`l_{II}`$ = 72.0, 78.0, and 84.0 degrees, respectively (Harmon, B.A. et al. 2002). Clearly interfering sources in various regions of the sky produce both positive and negative offsets in flux. These fluxes may be modulated at the precession period and may change sign, depending on the limb geometry.

Data gaps arise due to three specific causes: passage through the South Atlantic Anomaly when the detectors are switched off; known interfering source occulting within 10 seconds of the target; periods when the source is occulted by the earth’s limb at too shallow an angle for a significant detection to be made (or is not occulted at all). The latter effect happens only for sources lying at greater than 41 degrees in declination.

Interfering sources that are positively identified should be fitted for in the occultation technique, however this is (as in the case of the background correction) a question of optimisation and so interfering sources are not always completely removed. A signature of source interference is significant negative flux indicating that a persistent source close to the target has not been fitted correctly in calculating the occultation step-height for the target source. Provided the interfering source has not varied on timescales comparable to the recurrence period of the target, it will be averaged out in the epoch folding procedure. If however the source has varied then its effects must be positively identified before a folded lightcurve can be produced.

Unidentified soures of interference pose an additional problem in that their position relative to the earth’s limb and the target change continuously. This leads to modulation at the precession period as the interference contribution varies cyclically with the orbital plane.This type of artefact may be non sinusoidal although subtracting a sinewave with correct period and amplitude can attenuate it in some cases. In principle such signals could be used iteratively to improve the occultation analysis but this is not practical at this time. Persistent and identified transient sources of interference can be investigated on a source by source basis. In some cases datapoints or stretches of data have to be excluded, if for example a nearby source flares unpredictably. While in some datasets the interference is clearly apparent as a completely distinct feature in the power spectrum and can be attenuated by subtracting a sinewave fit to the modulation if neccesary. Scargle \[Scargle 1989\] states that the modified power spectrum (Lomb-Scargle periodogram) is equivalent to least-squares fitting with sinewaves, hence our procedure is equivalent to filtering the data in the frequency domain.

### 2.1 Implications of uneven sampling

Unevenly sampled data cannot be correctly analysed by the direct Fourier Transform technique as the resulting power spectrum is a convolution of both the source frequencies and sampling frequencies. This leads to spurious frequencies and possible suppression of genuine signals. Specifically; uneven sampling increases leakage of power to nearby frequencies, although at the same time, aliasing may be reduced. Our data are very unevenly sampled, especially when we excluded subsets of the data e.g. data intervals selected on the basis of an ephemeris in order to search for the possible 16.65d transient modulation in A0538-668. Fortunately the Lomb-Scargle periodogram is formulated for exactly this situation and is relatively insensitive to the sampling structure (Scargle 1982, 1989).

A further complication is the non-sinusoidal nature of many of the lightcurves, generally the less sinusoidal a signal is, the more components it will have in any frequency space representation. These components are seen in power spectra as sidebands to the main frequency spike. If the shape of an otherwise regular modulation changes with time, these peaks can become broadened, shift and become blended together. Thus in our period determinations we backed up the Lomb-Scargle result with Phase Dispersion Minimisation analysis (PDM). \[Stellingwerf 1978\]. Such periodograms typically contain a number of local minima corresponding to (For example) integer multiples of the true period. However the PDM technique is totally insensitive to the shape of the lightcurve and this makes for highly accurate period determinations in the case of irregular lightcurves. Periods were measured from the PDM periodogram by calculating the centroid of the feature of interest. We relied on the Lomb-Scargle periodogram for our primary detection tool because the separation of Fourier components is better suited to identification of interfering sources, sidebands, aliasing and statistical modeling of the signal. PDM was implemented in the Starlink program ’Period’.

### 2.2 Evaluation of signal detection confidence levels

In detecting weak, unevenly sampled, periodic signals against a varying background, statistical significance is a key issue. The only reliable way to evaluate the confidence levels associated with a power spectrum under these conditions is by statistical simulation.

We first measured the mean and variance of each dataset, we then created a series of random numbers conforming to a Gaussian distribution with the same mean and variance as each dataset. This series was then combined with the sampling structure (Timing information) of the original dataset by simply replacing each flux measurement with a number drawn from the random gaussian distribution. Thus creating a simulation of what BATSE might have seen, had there been no periodic X-ray sources in its field of view. This dataset was then subjected to Lomb-Scargle analysis and the power of the highest peak and the highest peak in a specified frequency range recorded. This process was then repeated 1000 times after which the mean and and variance of the peak noise power (PNP)were calculated, both for the whole spectral range available and for the specified frequency range of interest. For the Lomb-Scargle periodograms presented here, the displayed confidence levels are the mean PNP + 3 standard deviations. Thus our confidence levels directly reflect how likely it is that a given signal is real or due to chance without making any assumptions about the underlying properties of the data. The confidence levels calculated by this procedure for each dataset were later used in evaluating the Lomb-Scargle periodograms generated from simulated datasets.

### 2.3 Jitter in outburst turn-on times

The width (FWHM) of peaks in the power spectrum may be related to orbit-to-orbit variation in the exact time at which outbursts occur. We investigated this by comparing the FWHM of the orbital modulation peak in the Lomb-Scargle periodogram to that measured for a simulated dataset containing a sinewave of appropriate period and amplitude superimposed on a background created by randomizing the flux readings in the original dataset. Results are shown in Table 1 for comparison. Further investigation of this parameter by a Monte-Carlo technique showed that the FWHM is dependent on the noise present. Specifically the FWHM becomes more variable as the S/N ratio deteriorates. Combined with the well known dependence of the FWHM on the length of the dataset, the resolution of the periodogram in respect of this parameter had to be calculated independently for each dataset. For each source a Monte-Carlo simulation was run by first finding the correct amplitude sinewave to give the same S/N ratio as that seen in the data -using the limit finding procedure described in section 3.1.1. The flux values in the data were then repeatedly randomized and added to this signal. The Lomb-Scargle periodogram was calculated each time and the FWHM of the peak measured and recorded. After 100 cycles the mean and standard deviation of the FWHM were calculated and the results are displayed in Table 1.

In addition to movement in the outburst phase it is possible that other effects such as variations in the shape of the outburst profile might also give rise to a widening of the peak in the power spectrum. However, inspection of the results presented in Table 1 shows no evidence of any variability in the power spectrum peak indicating a general lack of jitter in turn-on times and no significant variation in the outburst profiles.

## 3 Results

Table 1 summarises key observational parameters of the HMXBs studied, while Figure 1 shows the X-ray lightcurves folded over 20 phase-bins. The errors shown are observational errors inherent in the occultation analysis, in producing the averaged lightcurves, the error-bars are the standard deviation of points in each phase-bin.

### 3.1  Link between Type I X-ray outbursts and long period optical variations in A0538-668

We tested the prediction that A0538-668 only exhibits X-ray outbursts correlated with its orbital period at times coincident with minima in its long period optical lightcurve. \[Alcock et al 1999\] We note that previous widely spaced X-ray observations of A0538-668 tend to support this hypothesis, however we find conclusively that between TJD 8363 - 10999 (The time span of the BATSE occultation data) no significant modulation of 20-70 keV X-ray flux at the orbital period is present at the predicted epochs. To do this we calculated an ephemeris for precise times at which this model predicts outbursts to occur. The orbital parameters were taken from Skinner \[Skinner et al 1982\], P=16.65d, with a zeropoint T0=TJD 3423.96 Maxima are expected at phase 0.9 in this scheme if the X-rays follow the optical variability. The long period parameters from Alcock et al \[Alcock et al 1999\] are under this scheme: P=420.82d, with arbitrary zeropoint TJD 9001. The long period minima occur at arbitrary phase 0.55. The ephemeris was then used to isolate time intervals when X-ray outbursts were expected. The duration of the time interval examined either side of long period minima was 100 days to cover the data coincident with the optical lightcurve ”dip” regions.

According to Skinner \[Skinner et al 1982\], A0538-668 was in the field of view of the Einstein satellite on four occasions. According to the ephemeris described here, detections reported on 16 Dec 1980 and 3 Feb 1981 lie in the optical dip, and non detections on 10 Apr 1979 and 30 May 1980 lie firmly in the centre of the bright phase of the optical lightcurve.

#### 3.1.1 BATSE detection threshold for A0538-668

By creating test datasets with the same mean and variance as the BATSE data and containing only simulated signals, we determined how strong a periodic signal must be in order to give a 3$`\sigma `$ detection. This approach could also be used to model other components of a power spectrum, particularly where interfering sources and sampling structure lead to spurious frequencies.

Known interfering signals in the A0538-668 data include LMC X-4 (which is only 0.6$`\mathrm{°}`$ away and dominates the periodic flux in this dataset) and the spacecraft precessional modulation. These components were modelled by measuring their amplitudes in folded lightcurves at 35.5d and 52d respectively.

Using the ephemeris described above BATSE data were selected which were coincident with the four long-period optical minima observed during the MACHO observations of Alcock et al \[Alcock et al 1999\]. This resulted in four blocks of BATSE one-day averages each of duration 100 days spanning the time interval TJD 9124.52 - TJD 10999.5 Having failed to detect any statistically significant modulation consistent with 16.65 days, this subset of the BATSE observations of A0538-668 was then repeatedly randomized to destroy any periodic content and combined with an incrementaly increasing sinusoidal signal plus fixed sinewave components to simulate the interfering sources in the real data. After each step, the Lomb-Scargle periodogram was computed and the highest peak in a pre-determined frequency range was compared to the 3 $`\sigma `$ confidence level previously calculated for that frequency range as described in section 2.3. Once this threshold was exceeded, the amplitude of the sine wave signal was recorded. The detectability of such a weak signal is somewhat influenced by the noise present. At the very margins of detectability the power in a given frequency can fluctuate up and down despite a gradual increase in signal amplitude. Hence we used the criterion that the power in the selected frequency range must remain above the threshold value irrespective of chance structures in noise in order to qualify. The derived upper limit was independently tested using many different noise sets to verify this criterion was met, i.e. That a signal of the amplitude quoted would definitely have been detected by the methods described here, had it been present. Weaker signals would not consistently rise above the threshold although detections of marginal significance are possible.

## 4 Discussion

### 4.1 EXO2030+375

This object has been subject of several recent studies. Reig et al \[Reig et al 1998\] presented results of multi-wavelength observations and modeling of the accretion rate. Wilson et al (2002) report on a decade of X-ray observations combined with IR and H$`\alpha `$ measurements which show a decline in the density of the circumstellar disk around the Be star. This decline was followed by a sudden drop in the X-ray flux and a turn-over from a spin-up trend to spin-down in the frequency history. This is the first Be/X-ray binary to show an extended interval, about 2.5 years, where the global trend is spin-down, but the outbursts continue. In 1995 the orbital phase of EXO 2030+375’s outbursts shifted from peaking about 6 days after periastron to peaking before periastron. The outburst phase slowly recovered to peaking at about 2.5 days after periastron. Wilson et al (2002) interpret this shift in orbital phase followed by a slow recovery as evidence for a global one-armed oscillation propagating in the Be disk. This was further supported by changes in the shape of the H$`\alpha `$ profile which are commonly believed to be produced by a reconfiguration of the Be disk. The truncated viscous decretion disk model provided an explanation for the long series of normal outbursts and the evidence for an accretion disk in the brighter normal outbursts.

In the present analysis we find convincing evidence for elevated flux levels at apastron in the earlier epoch (see Figure 1). The RXTE ASM lightcurve for epoch 2 was studied by Reig et al \[Reig et al 1998\], we also analysed this data in order to compare the lightcurve and power spectrum with the BATSE results. The power spectra from both instruments appear to share the same features in epochs 1 & 2, peaks at $`P_{orb}`$ & $`P_{orb}/2`$ with powers at similar ratios. Direct comparison of the BATSE and RXTE folded lightcurves suggests that X-ray emission from EXO2030+375 may be softer during the apastron outbursts than during the normal Type I outbursts at periastron. This hypothesis is drawn from the relative amplitudes of the two peaks in the lightcurve and the differing spectral responses of the two instruments.

Comparison of the FWHM of the orbital period peak at 46d with that derived from simulated data reveals that although the phase of peak outburst intensity shifts between the two epochs, its recurrence is otherwise periodic. There is no evidence for any gradual shift in phase in either epoch.

### 4.2 A0538-668

We find that during the period coincident with the MACHO monitoring of Alcock et al \[Alcock et al 1999\], A0538-668 did not exhibit hard X-ray outbursts at the 16.65 day orbital period during either the optical minima or maxima. The upper limit for such activity is found to be 3.0 $`\times 10^3`$ photon cm-2 s-1at 20-70 kev. We compare this limit with the X-ray spectrum obtained with Ariel V \[White & Carpenter 1978\] at the peak of the July 1977 outburst. By assuming a mean photon energy of 45 keV over a spectral range of 20-70 keV for the BATSE observations we estimate a maximum flux density of 2.7 $`\times 10^3`$ keV cm-2 s-1 keV-1. The fitted spectrum of White & Carpenter \[White & Carpenter 1978\] is extremely steep and dominated by a thermal component. By extrapolating their model, our lower limit for detection by BATSE in the data presented here is consistent with the flux that would be expected on the basis of this model to within errors. Thus we conclude that had A0538-668 been active at its previous level BATSE would have detected it.

### 4.3 Centaurus X-3

This source was detected in BATSE data at 29$`\sigma `$ with a period 2.087d. Cen X-3 is the only HMXB to have a well determined orbital period derivative \[Nagase et al 1992\] however this is too small to show up in our analysis. The power spectrum in Figure 2 shows two peaks, but comparison with the power spectrum for simulated data showed the higher-frequency peak to be an artefact of aliasing. The folded lightcurve Figure 1 is a simple sawtooth like form, brightening more rapidly than it decays.

### 4.4 4U1907+097

Two peaks with $``$3$`\sigma `$ significance are found at periods of 16.8d and 7.39d. The orbital period of the system is believed to be 8.38d (Marshall and Ricketts 1980, Cook and Page 1987) thus the 16.8d X-ray modulation is consistent with twice the binary period. However, there is no evidence for a peak in the power spectrum at $``$8.4d - this is not surprising if one considers the folded lightcurve in Figure 1 which shows just a single broad peak when folded at 16.8d. Constructing a simulated dataset containing a sinewave at 8.38d with amplitude taken as 0.002 (measured from a lightcurve folded at this period) failed to reproduce the 16.8d modulation. Thus it appears that the emission is really modulated at twice the binary period in these BATSE data.

### 4.5 4U1700-377

Early BATSE occulation data from this source have already been the subject of a paper by Rubin et al., 1996. Our PDM period agrees with 3.41180$`\pm `$0.00003d from Copernicus X-ray observation \[Braduardi 1978\]. We also detect significant peaks in the power spectrum at 13.808d at a significance of 11.3$`\sigma `$, and at exactly half this period. This period has been previously reported by Konig & Maisack \[Konig & Maisack 1997\]. We further investigated this period by subtracting the orbital modulation from the lightcurve, using Figure 1 as a template. This procedure resulted in reduction of spectral power at the orbital period by a factor of 24. No attenuation of the 13.808d period was observed, in fact a significant increase resulted. Hence we are confident that the 13.808d period is an independent modulation and not an harmonic of the orbital period.

### 4.6 4U1145-619

A peak was detected in LS periodogram at 190.238d with FWHM=17.48 days. A PDM periodogram determination of 186.6d $`\pm `$0.348068d measured by centroiding is consistent with the established value of 186.5 days.

### 4.7 OAO1657-41

The PDM period of 10.4626$`\pm `$0.0011d is in agreement with the period of 10.4436$`\pm `$0.0038d \[Chakrabarty et al 1993\] from BATSE pulse timing studies.

### 4.8 LMC X-4

The 1.4d binary period is not detected. The ”average” Nyquist period of the daily average data is 2.4d so the single-step data was also analysed, again with no detection. The 30.5d X-ray modulation believed to be due to the precession of a warped accretion disk \[Lang et al 1981\] is detected and is shown to be persistent throughout the full span of the data Figure 1. All period finding techniques tried, (Lomb-Scargle, PDM, CLEAN) derive a period for this modulation slightly shorter than usually quoted. A PDM periodogram constructed at a resolution of 2$`\times 10^5`$ cycles and measured by centroiding, derived a period of 30.349d while the L-S periodogram yielded 30.355d. Monte Carlo simulations of sinusoidal signals at 30.5d and 30.35d superimposed on noise created by randomizing the original dataset show that the LS periodogram frequency resolution is 1.43765$`\times 10^5`$ at the S/N level seen in the data. A solution at 30.5d is thus seemingly in conflict with the BATSE dataset. A remarkable fact is the stability of this modulation over the 7 years of data, since no jitter is evident from the values in Table 1

### 4.9 GX301-2

This system is clearly the brightest, most regular source studied here. The phase-binned folded lightcurve shows increasing flux levels consistent with phase 0.5, each of which reached type-I outburst levels. This result confirms previous work by Pravdo et al, 1995 and Koh et al, 1997 - both of these works also present evidence for apastron outbursts from this system. The power spectrum for GX301-2 shows a pattern of strong peaks which after comparison to simulated data are found to be real. These frequency components, appearing at integer submultiples of the orbital period are a signature of the highly non-sinusoidal lightcurve. The spectrum shown is that generated after the removal of 4 particularly bright apastron flux points, as expected their inclusion boosts the $`P_{orb}/2`$ component.

## 5 Conclusions

Observations of accreting pulsars in high mass binary systems using the BATSE occultation technique reveal that this approach is sufficiently sensitive to derive accurate measurements of the X-ray modulation. Evidence for apastron X-ray outbursts is found in 2 sources: EXO2030+375 and 4U1145-619. All of these sources are in high mass, long period systems, EXO2030+375 and 4U1145-619 have similarities in their folded lightcurves having flux levels significantly above minimum for most of each orbit and apastron features reaching 30-40 percent of peak flux. GX301-2 has a much flatter minimum level between emissions and the apastron emission is only 10 percent of peak flux. These differences undoubtedly reflect the different nature of the companion star, whether Be or supergiant. The fact that none of the shorter period systems show evidence for this phenomenon may be explainable in terms of the standard Bondi-Hoyle accretion model if suitable parameters can be estimated/calculated for the outflow velocity and density of the circumstellar material/wind along the orbital trajectory of the neutron star. We hypothesise that in short period systems the radial profile of the stellar wind can never satisfy the required conditions for apastron outbursts due to the smaller apastron separation distances and higher neutron star velocities involved. The FWHM of the orbital period peaks in the power spectra were investigated as a possible indicator of variability in the exact time at which outbursts ”switch on”. To within the uncertainties imposed by the signal/noise ratio of the data this effect was not detected in any of the sources.

## Acknowledgments

SL acknowledges the use of STARLINK software and the support of a PPARC studentship.


format: text

text:
# Modular functions and Ramanujan sums for the analysis of 1/𝑓 noise in electronic circuits

## I Introduction

$`1/f`$ noise, already discovered by Nyquist in resistors at the dawn of electronic age, still fascinates experimentalists and theoreticians. This is because none physical principle or mathematical function allows to predict easily the observed hyperbolic $`1/f`$ power spectral density. Nevertheless the occurrence of $`1/f`$ fluctuations in areas as diverse as electronics, chemistry, biology, cognition or geology claims for an unifying mathematical principle WentianLi .

A newly discovered clue for $`1/f`$ noise lies in the concept of a phase locked loop (or PLL) APL02 . In essence two interacting oscillators, whatever their origin, attempt to cooperate by locking their frequency and their phase. They can do it by exchanging continuously tiny amounts of energy, so that both the coupling coefficient and the beat frequency should fluctuate in time around their average value. Correlations between amplitude and frequency noise were observed Yamoto .

Fortunately one can gain a good level of understanding of phase locking from quartz crystal oscillators used in high frequency synthesizers, ultrastable clocks and communication engineering (mobile phones or so). The PLL used in a FM radio receiver is a genuine generator of $`1/f`$ noise. Close to phase locking the level of $`1/f`$ noise scales approximately as $`\stackrel{~}{\sigma }^2`$, where $`\stackrel{~}{\sigma }=\sigma K/\stackrel{~}{\omega _B}`$ is the ratio between the open loop gain $`K`$ and the beat frequency $`\stackrel{~}{\omega _B}`$ times a constant coefficient $`\sigma `$. The relation above is explained from a simple non linear model of the PLL known as Adler’s equation

$$\dot{\theta }(t)+KH(P)\mathrm{sin}\theta (t)=\omega _B,$$
(1)
where at this stage $`H(P)=1`$, $`\omega _B=\omega (t)\omega _0`$ is the angular frequency shift between the two quartz oscillators at the input of the non linear device (a Schottky diodes double balanced mixer), and $`\theta (t)`$ is the phase shift of the two oscillators versus time t. Solving (1) and differentiating one gets the observed noise level $`\stackrel{~}{\sigma }`$ versus the one $`\sigma =\delta \omega _B/\stackrel{~}{\omega _B}`$ for the open loop case. Thus the model doesn’t explain the existence of $`1/f`$ noise but correctly predicts its dependance on the physical parameters of the loop APL02 .

## II Low Pass Filtering and $`1/f`$ Noise

Besides one can get detailed knowledge of harmonic conversions in the PLL by accounting for the transfer function $`H(P)`$, where $`P=\frac{d}{dt}`$ is the Laplace operator. If $`H(P)`$ is a low pass filtering function with cut-off frequency $`f_c`$, the frequency at the output of the mixer + filter stage is such that

$$\mu =f_B(t)/f_0=q_i|\nu p_i/q_i|f_c/f_0,p_i\text{ and }q_i\text{ integers}.$$
(2)
The beat frequency $`f_B(t)`$ results from the continued fraction expansion of the input frequency ratio

$$\nu =f(t)/f_0=[a_0;a_1,a_2,\mathrm{}a_i,a,\mathrm{}],$$
(3)
where the brackets mean expansions $`a_0+1/(a_1+1/(a_2+1/\mathrm{}+1/(a_i+1/(a+\mathrm{}))))`$. The truncation at the integer part $`a=[\frac{f_0}{f_cq_i}]`$ defines the edges of the basin; they are located at $`\nu _1=[a_0;a_1,a_2,\mathrm{},a_i,a]`$ and $`\nu _2=[a_0;a_1,a_2,\mathrm{}a_i1,1,a]`$ Fluc01 . The two expansions in $`\nu _1`$ and $`\nu _2`$, prior to the last filtering partial quotient $`a`$, are the two allowed ones for a rational number. The convergents $`p_i/q_i`$ at level $`i`$ are obtained using the matrix products

$$\left[\begin{array}{cc}a_0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{cc}a_1& 1\\ 1& 0\end{array}\right]\mathrm{}\left[\begin{array}{cc}a_i& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}p_i& p_{i1}\\ q_i& q_{i1}\end{array}\right].$$
(4)
Using (4) one can get the fractions $`\nu _1`$ and $`\nu _2`$ as $`\nu _1=\frac{p_a}{q_a}`$ and $`\nu _2=\frac{p_i(2a+1)p_a}{q_i(2a+1)q_a}`$ so that, with the relation relating convergents $`(p_iq_{i1}p_{i1}q_i)=(1)^{i1}`$, the width of the basin of index $`i`$ is $`|\nu _1\nu _2|=\frac{2a+1}{q_a(q_a+(2a+1)q_i)}\frac{1}{q_aq_i}`$ whenever $`a>1`$ (see also Fig. 1).

In our previous publications we proposed a phenomenological model for $`1/f`$ noise in the PLL, based on an arithmetical function which is a logarithmic coding for prime numbers APL02 ,Fluc01 . If one accepts a coupling coefficient evolving discontinuously versus the time $`n`$ as $`K=K_0\mathrm{\Lambda }(n)`$, with $`\mathrm{\Lambda }(n)`$ the Mangoldt function which is $`\mathrm{ln}(p)`$ if $`n`$ is the power of a prime number $`p`$ and $`0`$ otherwise, then the average coupling coefficient is $`K_0`$ and there is an arithmetical fluctuation $`ϵ(t)`$

$`\psi (t)={\displaystyle \underset{n=1}{\overset{t}{}}}\mathrm{\Lambda }(n)=t(1+ϵ(t)),`$
$`tϵ(t)=\mathrm{ln}(2\pi ){\displaystyle \frac{1}{2}}\mathrm{ln}(1t^2){\displaystyle \underset{\rho }{}}{\displaystyle \frac{t^\rho }{\rho }}.`$ (5)
The three terms at the right hand side of $`tϵ(t)`$ come from the singularities of the Riemann zeta function $`\zeta (s)`$, that are the pole at $`s=1`$, the trivial zeros at $`s=2l`$, $`l`$ integer, and the zeros on the critical line $`\mathrm{}(s)=\frac{1}{2}`$Fluc01 . Also the power spectral density roughly shows a $`1/f`$ dependance versus the Fourier frequency $`f`$. This is an unexpected relation between Riemann zeros (the unsolved Riemann hypothesis is that all zeros should lie on the critical line) and $`1/f`$ noise.

We improved the model by replacing the Mangoldt function by its modified form $`b(n)=\mathrm{\Lambda }(n)\varphi (n)/n`$, with $`\varphi (n)`$ the Euler (totient) function PRE02 . This seemingly insignificant change was introduced by Hardy Hardy  in the context of Ramanujan sums for the Goldbach conjecture and resurrected recently by Gadiyar and Padma for analyzing the distribution of prime pairs Gadiyar . Then by defining the error term $`ϵ_B(t)`$ from

$$B(t)=\underset{n=1}{\overset{t}{}}b(n)=t(1+ϵ_B(t)),$$
(6)
its power spectral density $`S_B(f)\frac{1}{f^{2\alpha }}`$ exhibits a slope close to twice the Golden ratio $`\alpha (\sqrt{5}1)/20.618`$ (see Fig. 2).

The modified Mangoldt function occurs in a natural way from the logarithmic derivative of the quotient

$$Z(s)=\frac{\zeta (s)}{\zeta (s+1)}=\underset{n1}{}\frac{\varphi (n)}{n^{s+1}},$$
(7)
since $`\frac{Z^{}(s)}{Z(s)}=_{n1}\frac{b(n)}{n^s}`$. This replaces the relation from the Riemann zeta function where $`\frac{\zeta ^{}(s)}{\zeta (s)}=_{n1}\frac{\mathrm{\Lambda }(n)}{n^s}`$.

## III The Hyperbolic Geometry of Phase Noise and $`1/f`$ Frequency Noise

The whole theory can be justified by studying the noise in the half plane $`H=\{z=\nu +\mathrm{𝐼𝑦},I^\mathit{2}=\mathit{1},y>\mathit{0}\}`$ of coordinates $`\nu =\frac{f}{f_0}`$ and $`y=\frac{f_B}{f_c}>0`$ and by introducing the modular transformations

$$z\gamma (z)=\frac{p_iz+p_i^{}}{q_iz+q_i^{}},p_iq_i^{}p_i^{}q_i=1.$$
(8)
The set of images of the filtering line $`z=\nu +I`$ under all modular transformations can be written as

$$|z(\frac{p_i}{q_i}+\frac{I}{2q_i^2})|=\frac{1}{2q_i^2}.$$
(9)
Equation (9) defines Ford circles Rademacher Planat2 (see Fig. 3) centered at points $`z=\frac{p_i}{q_i}+\frac{I}{2q_i^2}`$ with radius $`\frac{1}{2q_i^2}`$ Rademacher . To each $`\frac{p_i}{q_i}`$ belongs a Ford circle in the upper half plane, which is tangent to the real axis at $`\nu =\frac{p_i}{q_i}`$. Ford circles never intersect: they are tangent to each other if and only if they belong to fractions which are adjacent in the Farey sequence $`\frac{0}{1}<\mathrm{}\frac{p_1}{q_1}<\frac{p_1+p_2}{q_1+q_2}<\frac{q_2}{q_2}\mathrm{}<\frac{1}{1}`$.

For modular transformations one easily calculates $`\mathrm{}\gamma (z)=\frac{y}{|q_iz+q_i^{}|^2},\frac{d\gamma (z)}{dz}=\frac{1}{(q_iz+q_i^{})^2}`$. where the symbol $`\mathrm{}`$ means the imaginary part. This implies $`|\frac{d\gamma (z)}{dz}|(\mathrm{}\gamma (z))^1=(\mathrm{}z)^1`$ and the invariance of the non-Euclidean metric $`\frac{|dz|}{y}=\frac{(d\nu ^2+dy^2)^{1/2}}{y}`$.

A basic fact about the modular transformations (8) is that they form a discontinuous group $`\mathrm{\Gamma }SL(2,Z)/\{\pm \mathit{1}\}`$, which is called the modular group. The action of $`\mathrm{\Gamma }`$ on the half-plane $`H`$ looks like the one generated by two independent linear translations on the Euclidean plane, which is equivalent to tesselate the complex plane $`C`$ with congruent parallelograms. Here one introduces the fundamental domain of $`\mathrm{\Gamma }`$ (or modular surface) $`F=\{zH:|z|\mathit{1},|\nu |\frac{\mathit{1}}{\mathit{2}}\}`$, and the family of domains $`\{\gamma (F),\gamma \Gamma \}`$ induces a tesselation of $`H`$Serre ,Gutzwiller .

Corresponding to the Riemannian metric there is the non-Euclidean Laplacian $`\mathrm{\Delta }=y^2(\frac{^2}{\nu ^2}+\frac{^2}{y^2})`$, and one can write a non-Euclidean wave equation Gutzwiller

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}(\mathrm{\Delta }+\frac{1}{4})\mathrm{\Psi }.$$
(10)
Stationary solutions $`\mathrm{\Psi }=\varphi (\nu ,y)\mathrm{exp}(I\omega t)`$ satisfy

$$(\mathrm{\Delta }+\frac{1}{4}+k^2)\varphi =0,\text{with}k^2=\frac{2mE}{\mathrm{}^2},E=\mathrm{}\omega .$$
(11)
They are special $`\nu `$-independent solutions $`\varphi =y^s\text{with}s=\frac{1}{2}+Ik`$, and any other wave can be transformed by (8) to form the new solution

$$\varphi _s(z)=(\mathrm{}\gamma (z))^s=\frac{y^s}{|q_iz+q_i^{}|^{2s}},s=\frac{1}{2}+Ik.$$
(12)
In Gutzwiller’s paper Gutzwiller  there is a detailed discussion on the geometry of solutions (12). The $`\nu `$-independent solution represents a wave propagating in the y-direction $`\mathrm{\Psi }(y)=y^{1/2}\mathrm{exp}(Ik\mathrm{ln}yI\omega t)`$. The factor $`y^{1/2}`$ comes from restricting the total flux in a vertical strip of constant width such as $`0<\nu <1`$, so that with the hyperbolic metric one has $`_0^1𝑑\nu /y=1/y`$. The phase factor is explained by looking at wavefronts $`y=C^{\mathrm{te}}`$ and the distance proportional to them is $`_1^y𝑑y/y=\mathrm{ln}y`$.

For equation (12) the wavefronts are obtained by assigning to the quantity $`\frac{y}{|q_iz+q_i^{}|^2}`$ some constant value that one can choose to be the unity. These wavefronts are circles of equation $`(\nu +\frac{q_i^{}}{q_i})^2+(y\frac{1}{2q_i^2})^2=\frac{1}{4q_i^4}`$. They are tangent to the real axis at $`\nu _0=\frac{q_i^{}}{q_i}`$ and of radius $`\frac{1}{2q_i^2}`$. Up to a constant shift along the real axis $`\nu `$, they are the same as the Ford circles shown on Fig. 3.

Thus an outgoing plane wave can be said to start at $`y=1`$, and move up into $`H`$ with increasing $`y`$, its wavefronts parallel to the real axis. All other plane waves start at some circle which touches the real axis at some arbitrary point $`\nu _0`$ of radius $`\frac{1}{2q_i^2}`$. Then they contract by shrinking their radius while maintaining their point of contact $`\nu _0`$. These wave fronts are horocycles in the hyperbolic geometry of the half-plane. They are perpendicular to the geodesics, which are Euclidean half-circles, hitting the real axis at $`\nu _0`$, at a right angle.

In addition to (12) they are general solutions Terras  of the form

$$\varphi _s(z)=\underset{\gamma \mathrm{\Gamma }_{\mathrm{}}/\mathrm{\Gamma }}{}(\mathrm{}\gamma (z))^s,s=\frac{1}{2}+Ik,$$
(13)
where $`\mathrm{\Gamma }_{\mathrm{}}`$ is the stabilizer in the modular group $`\mathrm{\Gamma }`$ of the cusp at infinity, i.e. the subgroup of integer translations $`zz+n`$. The solution corresponds to waves scattered again and again from the modular surface $`F`$. For arbitrary $`s`$ the series (13) can be rewritten as $`\varphi _s(z)=\frac{1}{2}y^s_{(p_i,q_i)=1}\frac{1}{|p_iz+q_i^{}|^{2s}}`$ Terras , with the summation performed over all coprime numbers $`(p_i,q_i)=1`$. The convergence is ensured for $`\mathrm{}(s)>1`$.

The series (13) satisfies a functional equation

$$\xi (2s)\varphi _s(z)=\xi (22s)\varphi _{1s}(z),$$
(14)
where $`\xi (s)=\pi ^{s/2}\mathrm{\Gamma }(s/2)\zeta (s)`$ is the completed Riemann zeta function. It follows that (13) is expanded as

$`\varphi _s(z)=y^s+S(s)y^{1s}+T_s(y),`$
$`\text{with}S(s)=A(s)Z(s),`$
$`Z(s)={\displaystyle \frac{\zeta (s1)}{\zeta (s)}},A(s)={\displaystyle \frac{\mathrm{\Gamma }(1/2)\mathrm{\Gamma }(s1/2)}{\mathrm{\Gamma }(s)}},`$ (15)
and the remaining term vanishes exponentially when $`y\mathrm{}`$. For $`s=\frac{1}{2}+Ik`$, that is $`s`$ on the critical line, the remainder is $`T_s(y)=0`$. In such a case the scattering coefficient equals

$$S(\frac{1}{2}+Ik)=\frac{\xi (2Ik)}{\xi (1+2Ik)},\text{with}|S(\frac{1}{2}+Ik)|=1,$$
(16)
that is the flux of the reflected wave is equal and opposite sign the incoming flux. The phase angle represented on Fig. 4 is a very complicated function of $`k`$ and is even considered as a prototype for quantum chaos Gutzwiller .

It is observed that the term $`A(\frac{1}{2}+Ik)`$ in (15) increases while $`|Z(\frac{1}{2}+Ik)|`$ decreases monotonously to preserve the modulus 1 of $`S(\frac{1}{2}+Ik)`$. All the stochastic dynamics is encoded in the function $`Z(\frac{1}{2}+Ik)`$. Fig. 4 represents the phase angle $`\theta (k)`$ in the scattering coefficient $`S(\frac{1}{2}+Ik)=\mathrm{exp}(I\theta (k))`$ in comparison to the one $`\theta _0(k)`$ of $`Z(\frac{1}{2}+Ik)`$. Apart for weak changes the two curves looks similar; this is confirmed from the power spectral density which is about the same for the two cases.

Looking at the logarithmic derivative of the scattering coefficient $`S(s)=\mathrm{exp}(i\theta (s))`$, one can get the counting function

$$\theta ^{}(k)=\frac{d\mathrm{ln}S(s)}{ds}\text{at}s=\frac{1}{2}+Ik.$$
(17)
which is related to the stochastic factor $`Z(s)`$ as

$$\frac{Z^{}(s)}{Z(s)}=\underset{n=1}{\overset{\mathrm{}}{}}\frac{b(n)}{n^s}=s_1^{\mathrm{}}t^{s1}B(t)𝑑t,$$
(18)
with $`B(t)=_{n=1}^{\mathrm{}}b(n)`$, $`b(n)=\mathrm{\Lambda }(n)\varphi (n)/n`$. By inverting the Mellin transform in (15) one gets

$$B(t)=\frac{1}{2I\pi }_{\mathrm{}(s)I\mathrm{}}^{\mathrm{}(s)+I\mathrm{}}\frac{Z^{}(s)}{Z(s)}.\frac{t^s}{s}ds.$$
(19)
This extends the calculation performed in (5) for getting the error term in the summatory Mangoldt function $`\psi (t)`$. There the Riemann zeta function $`\zeta (s)`$ replaces the quotient $`Z(s)=\zeta (s)/\zeta (s+1)`$ given in (7). As reminded in Fig. 2, the error term in the summatory modified Mangoldt function is very close to a $`1/f^{2\alpha }`$ noise, with $`\alpha `$ the Golden ratio.

## IV Concluding Remarks

From its definition $`1/f`$ noise is attached to the use of the fast Fourier transform (FFT). But the FFT refers to the fast calculation of the discrete Fourier transform (DFT) with a finite period $`q=2^l`$, $`l`$ a positive integer. In the DFT one starts with all $`q^{\mathrm{th}}`$ roots of the unity $`\mathrm{exp}(2i\pi p/q)`$, $`p=1\mathrm{}q`$ and the signal analysis of the arithmetical sequence $`x(n)`$ is performed by projecting onto the $`n^{\mathrm{th}}`$ powers (or characters of Z/qZ) with well known formulas.

The signal analysis based on the DFT is not well suited to aperiodic sequences with many resonances (by nature a resonance is a primitive root of the unity: $`(p,q)=1`$), and the FFT may fail to discover the underlying structure in the spectrum. We recently introduced a new method based on Ramanujan sums PRE02 ,Hardy .

$$c_q(n)=\underset{p=1}{\overset{q}{}}\mathrm{exp}(2i\pi \frac{p}{q}n)\text{with}(p,q)=1,$$
(20)
which are $`n^{\mathrm{th}}`$ powers of the $`q^{\mathrm{th}}`$ primitive roots of the unity. The sums are evaluated from the use of Möbius transforms as $`c_q(n)=\mu \left(\frac{q}{(q,n)}\right)\frac{\varphi (q)}{\varphi \left(\frac{q}{(q,n)}\right)}`$. Here $`(q,n)`$ is the greatest common divisor $`(q,n)`$ of $`q`$ and $`n`$ and $`\mu (n)`$ is the Möbius function. It is $`0`$ if $`n`$ contains a square, $`1`$ if $`n=1`$ and $`(1)^k`$ if $`n`$ is the product of $`k`$ distinct primes.

The sums are quasiperiodic versus the time $`n`$ ($`c_1=\overline{1};c_2=\overline{1,1};c_3=\overline{1,1,2}`$, where the bar indicates the period; for example $`c_3(4)=1`$) and aperiodic versus the order $`q`$ of the resonance. In particular $`c_q(n)=\mu (q)`$ whenever $`(n,q)=1`$.

Möbius function can be considered as a coding sequence for prime numbers, as it is the case of Mangoldt function. Mangoldt function is related to Möbius function thanks to the Ramanujan sums expansion found by Hardy Gadiyar

$$b(n)=\frac{\varphi (n)}{n}\mathrm{\Lambda }(n)=\underset{q=1}{\overset{\mathrm{}}{}}\frac{\mu (q)}{\varphi (q)}c_q(n).$$
(21)
We call such a type of Fourier expansion a Ramanujan-Fourier transform (RFT). General formulas are given in our recent publication PRE02  and in the paper by Gadiyar Gadiyar . This author also reports on a stimulating conjecture relating the autocorrelation function of $`b(n)`$ and the problem of prime pairs. In the special case (21), it is clear that $`\mu (q)/\varphi (q)`$ is the RFT of the modified Mangoldt sequence $`b(n)`$.

Using Ramanujan-Fourier analysis the $`1/f^{2\alpha }`$ power spectrum gets replaced by a new signature shown on Fig. 5, very close to $`\mu (q)/\varphi (q)`$ (up to a scaling factor). There is thus a deep relationship between $`1/f`$ noise in phase locking, the Golden ratio, the Möbius function, the modified Mangoldt function, the frequency of windings around the modular surface, and the Riemann hypothesis. All these ingredients arise from the hyperbolic geometry of the half-plane.


text:
# Berry phase in Magnetic Superconductors

## Abstract

In magnetic systems, electronic bands often acquire nontrivial topological structure characterized by gauge flux distribution in momentum($`𝐤`$)-space. It sometimes follows that the phase of the wavefunctions cannot be defined uniquely over the whole Brillouin zone. In this Letter we develop a theory of superconductivity in the presence of this gauge flux both in two- and three-dimensional systems. It is found that the superconducting gap has “nodes” as a function of $`𝐤`$ where the Fermi surface is penetrated by a gauge string.

Geometric phases and topology of wavefunctions have been the subject of intensive studies berry ; review . The quantized Hall effect (QHE) in a 2D electron system under a strong magnetic field can be also interpreted in terms of the topological integer called Chern number tknn . If the system has a gap, the Hall conductivity is proportional to the Chern number and is quantized. This nontrivial topology of Bloch waves is not specific to the QHE, but is common and universal in systems with broken time-reversal symmetry. Namely, in magnetic materials with the spin-orbit interaction and/or tilting of spin configuration, there occur many band-crossing points which are locally described as Weyl fermions and act as (anti)monopoles in the momentum-space ohgushi ; onoda ; shindou ; taguchi . The existence of (anti)monopoles means that the Bloch wavefunction cannot be defined in a single gauge choice sakurai . The anomalous Hall effect in ferromagnets represents this topological property of the Bloch wavefunctions ohgushi ; taguchi .

On the other hand, coexistence of magnetic ordering and superconductivity (SC) is found in many materials and is a recent significant issue super . Although detailed analysis of each material is still missing in most of the cases, it is highly desirable to establish the general feature of the SC in systems with broken time-reversal symmetry.

In the present Letter we shall study the SC made out of the Bloch states with nontrivial topology. The quantity of central importance is the gauge flux defined in the momentum space, which is generated by the vector potential $`𝐀_{(n)}(𝐤)=i𝐤n|_𝐤|𝐤n`$, where $`|𝐤n`$ is the Bloch wavefunction of the $`n`$-th band. Band indices are represented by subcripts with parenthesis. This $`𝐀_{(n)}(𝐤)`$ represents an overlap of the two wavefunctions separated infinitestimally in the $`𝐤`$-space, i.e. the Berry phase connection. The gauge flux is defined as $`𝐁_{(n)}(𝐤)=_𝐤\times 𝐀_{(n)}(𝐤)`$. It is worth noting that the monopole density $`\rho _{(n)}(𝐤)=\frac{1}{2\pi }_𝐤𝐁_{(n)}(𝐤)`$ is nonzero, and is given by $`\rho _{(n)}(𝐤)=_lq_{ln}\delta (𝐤𝐤_{ln})`$, where $`q_{ln}`$ is an integer called the strength of the monopole berry . The monopoles are located at “diabolical points” berry ; volovik , where bands touch each other. The monopole and antimonopole act as the source and the sink of the gauge flux $`𝐁(𝐤)`$, respectively, as in Fig. 1(a). In general, the up- and down-spins have different band structures; therefore, $`𝐀`$, $`𝐁`$, and $`\rho `$ are spin-dependent, and we shall write as $`𝐀_\alpha `$, $`𝐁_\alpha `$ and $`\rho _\alpha `$, where $`\alpha `$ is a spin index.

The question we address below is an effect of this gauge flux $`𝐁(𝐤)`$ on the SC properties. In order to characterize topological properties of the gap function $`\mathrm{\Delta }(𝐤)`$, we define a vector potential $`𝐀^\mathrm{\Delta }(𝐤)=_𝐤\text{Im}\mathrm{ln}\mathrm{\Delta }(𝐤)=_𝐤\text{arg}\mathrm{\Delta }(𝐤)`$, a flux density $`𝐁^\mathrm{\Delta }=_𝐤\times 𝐀^\mathrm{\Delta }(𝐤)`$, and a monopole density $`\rho ^\mathrm{\Delta }(𝐤)=\frac{1}{2\pi }_𝐤𝐁^\mathrm{\Delta }(𝐤)`$. For the moment, we assume that $`\mathrm{\Delta }(𝐤)`$ is nonvanishing in almost all over the BZ noteDelta ; we shall discuss later what happens without this assumption. As explained later, distribution of $`𝐁^\mathrm{\Delta }(𝐤)`$ is confined into strings as in Fig. 1(b), due to bosonic nature of $`\mathrm{\Delta }(𝐤)`$. The main result of the present Letter is that in 3D topological structures of the gap function $`\mathrm{\Delta }(𝐤)`$ is solely determined from that of the wavefunction: $`\rho _{\alpha \beta }^\mathrm{\Delta }(𝐤)=\rho _\alpha (𝐤)\rho _\beta (𝐤)`$. It reflects the fact that $`\mathrm{\Delta }_{\alpha \beta }(𝐤)`$ represents a pairing between $`(𝐤,\alpha )`$ and $`(𝐤,\beta )`$ electrons. Thus if the band structure in the normal state is known, we can immediately calculate the monopole density $`\rho _{\alpha \beta }^\mathrm{\Delta }(𝐤)`$, irrespective of details of the attractive potential. The monopole density $`\rho _{\alpha \beta }^\mathrm{\Delta }(𝐤)`$ tells us about zeros and phases of the gap function.

Henceforth we focus on the SC where only one component of the gap function is nonzero: (i) singlet SC, (ii) triplet SC with $`\mathrm{\Delta }_{}0`$, $`\mathrm{\Delta }_{}=\mathrm{\Delta }_{}=0`$, and (iii) triplet SC with $`\mathrm{\Delta }_{}0`$, $`\mathrm{\Delta }_{}=\mathrm{\Delta }_{}=0`$ note-triplet . The case (iii) is appropriate for half-metallic magnetic SC. In the singlet (triplet) SC, we have $`\mathrm{\Delta }(𝐤)=\mathrm{\Delta }(𝐤)`$ ($`\mathrm{\Delta }(𝐤)=\mathrm{\Delta }(𝐤)`$), implying $`𝐁^\mathrm{\Delta }(𝐤)=𝐁^\mathrm{\Delta }(𝐤)`$, and $`\rho ^\mathrm{\Delta }(𝐤)=\rho ^\mathrm{\Delta }(𝐤)`$ note-delta . We note that the time-reversal symmetry breaking is a necessary condition for nontrivial topology of the gap function. In (iii), broken time-reversal symmetry is assumed from the outset. We have $`\rho _{}^\mathrm{\Delta }(𝐤)=\rho _{}(𝐤)\rho _{}(𝐤)`$. Thus, inversion symmetry must be broken to have nonzero $`\rho _{}^\mathrm{\Delta }`$. In (i) and in (ii) we have $`\rho _{}^\mathrm{\Delta }(𝐤)=\rho _{}(𝐤)\rho _{}(𝐤)`$. Hence in nonmagnetic SC, $`\rho _{}^\mathrm{\Delta }(𝐤)`$ vanishes. In magnetic SC, including both ferromagnets and antiferromagnets, $`\rho _{}^\mathrm{\Delta }(𝐤)`$ is nonzero.

Let us consider the distribution of $`𝐁^\mathrm{\Delta }(𝐤)`$. $`𝐁^\mathrm{\Delta }(𝐤)=_𝐤\times _𝐤\text{arg}\mathrm{\Delta }(𝐤)`$ is nonzero only when $`\mathrm{\Delta }=0`$, i.e. $`\text{Re}\mathrm{\Delta }=\text{Im}\mathrm{\Delta }=0`$, which determines a curve in the 3D BZ in general (Fig. 1 (b)). Distribution of $`𝐁^\mathrm{\Delta }(𝐤)`$ is like a Dirac string, starting from a monopole and terminating at an antimonopole. In contrast, $`𝐁(𝐤)`$ does not share this feature; distribution of $`𝐁(𝐤)`$ spreads over the BZ (Fig. 1 (a)). Now we discuss the relation among the Fermi surface (FS), the (anti)monopoles and the Dirac string. For any closed surface $`S`$ in the BZ, the total flux penetrating it is quantized. We take the FS as $`S`$ for gapless systems. When the total flux penetrating one of the FS’s is nonzero, the gap function cannot be expressed as continuous. It is quite contrary to our conventional view note-2 . Because $`\rho ^\mathrm{\Delta }(𝐤)=\rho ^\mathrm{\Delta }(𝐤)`$, a FS symmetric with respect to $`𝐤=0`$ encloses zero monopoles in total, and $`\mathrm{\Delta }(𝐤)`$ can be continuous on this FS (Fig. 2(a)). There are, however, two cases where the Dirac strings penetrate the FS’s. One is the case with broken inversion as well as time-reversal symmetry. Then, the FS is not symmetric with respect to $`𝐤=0`$, and can enclose monopoles with nonvanishing strength (Fig. 2 (b)). The other is a pair of FS’s symmetric with respect to $`𝐤=0`$. It is possible that one FS encompasses a monopole and the other encompasses an antimonopole (Fig. 2(c)). For example, if one of the FS encloses a monopole of unit strength, a Dirac string starting from this monopole necessarily penetrates the FS. An intersection of this curve with the FS is nothing but a point node. This string can intersect the FS more than once, and the number of point nodes on this FS should be odd. The nonvanishing total flux implies that $`\mathrm{\Delta }(𝐤)`$ cannot be a single continuous function on each FS. The FS should be divided into regions, in each of which $`\mathrm{\Delta }`$ is continuous.

Here we remark on “conventional” nodes in anisotropic SC. Line nodes appear when $`\text{Re}\mathrm{\Delta }`$ and $`\text{Im}\mathrm{\Delta }`$ have a common real factor $`f(𝐤)`$ as $`\mathrm{\Delta }(𝐤)=f(𝐤)[\text{Re}g(𝐤)+i\text{Im}g(𝐤)]`$. Zeros of $`f(𝐤)`$ determine a surface, whose intersection with the FS is a line node. It does not affect $`𝐁^\mathrm{\Delta }(𝐤)`$, because $`\text{Im}\mathrm{ln}\mathrm{\Delta }(𝐤)=\text{Im}\mathrm{ln}g(𝐤)`$. Hence, only point nodes are concluded from the gauge flux argument. Difference of nonmagnetic SC from magnetic SC lies only in absence of magnetic monopoles; in nonmagnetic SC, Dirac strings necessarily form loops, which may cross the boundary of the BZ. There are many experimental methods to get information on point nodes. In particular, heat conductivity measurements tell us about the direction of the nodes izawa . When the positions of the nodes are known, care must be taken to assign the phase of $`\mathrm{\Delta }(𝐤)`$ for magnetic SC, because our theory asserts that $`\mathrm{\Delta }(𝐤)`$ need not be continuous on the FS. Group-theoretical classifications of gap functions in su ; vg  do not apply to the SC with nontrivial topology, since they assume that gap functions are continuous. It would be interesting to classify all possible multivalued gap functions.

To see why topological structure of the wavefunction is inherited by the gap function, we note the following. If the wavefunction is a continuous function of $`𝐤`$, there is no monopole. Thus, if there is a monopole, we should divide a surface surrounding it into regions, in each of which the wavefunction is continuous kohmoto . This resembles Dirac monopoles wy , where the vector potential cannot be a single continuous function in the whole space. This “patch” structure of the wavefunctions remains in the gap function via the BCS term in the Hamiltonian.

To understand this mechanism in detail, we start with a 2D model and generalize it to 3D. Similar topological arguments as we have developed in 3D are also possible in 2D as well. In 2D, we define a Chern number and a total vorticity of $`\mathrm{\Delta }(𝐤)`$ as $`\text{Ch}=\frac{1}{2\pi }_{\text{BZ}}d^2𝐤B(𝐤)_z`$ and $`\nu =\frac{1}{2\pi }_{\text{BZ}}d^2𝐤B^\mathrm{\Delta }(𝐤)_z`$, respectively. We can show that the total vorticity $`\nu _{\alpha \beta }`$ of the gap function $`\mathrm{\Delta }_{\alpha \beta }`$ is the sum of the Chern numbers for spin $`\alpha `$ and for spin $`\beta `$: $`\nu _{\alpha \beta }=\text{Ch}_\alpha +\text{Ch}_\beta `$ note-2D . We shall explain how this result comes out in the singlet SC on the 2D model proposed by Haldane haldane  as an illustrative example. We take this model because it is the minimal model which exhibits nontrivial topological structure of wavefunctions without a uniform magnetic field. It is a tight-binding model on the honeycomb lattice. Because there are two sublattices, the Hamiltonian is written as a $`2\times 2`$ matrix, which is conveniently expressed with the Pauli matrices $`𝝈`$ as $`H(𝐤)=vI+𝐰𝝈`$, where $`v=2t_2\mathrm{cos}\varphi _j\mathrm{cos}(𝐤𝐛_j)`$, $`w_1=t_1_j\mathrm{cos}(𝐤𝐚_j)`$, $`w_2=t_1_j\mathrm{sin}(𝐤𝐚_j)`$, and $`w_3=M2t_2\mathrm{sin}\varphi _j\mathrm{sin}(𝐤𝐛_j)`$, where $`𝐚_1=a(1,0)`$, $`𝐚_2=a(\frac{1}{2},\frac{\sqrt{3}}{2})`$, $`𝐚_3=𝐚_1𝐚_2`$, $`𝐛_1=𝐚_2𝐚_3`$, $`𝐛_2=𝐚_3𝐚_1`$, $`𝐛_3=𝐚_1𝐚_2`$. $`t_1`$, $`t_2`$, $`M`$, $`\varphi `$ are constants. The Pauli matrices are not associated with spins, but with the sublattice structure. We set $`M=0`$ and $`0<\varphi <\pi `$. Time-reversal symmetry is broken by a staggered flux, while the uniform field is absent. This staggered flux is not necessarily external; staggered flux can be realized by a spontaneous magnetic moment in the material itself, through spin-chirality and/or spin-orbit coupling ohgushi ; taguchi . Let us calculate the Chern number for the lower band. By dividing the BZ into two regions $`V_\pm `$ as in Fig. 3, where $`𝐤_\pm ^0=\pm \frac{2\pi }{a}(\frac{1}{3},\frac{1}{3\sqrt{3}})`$ are two inequivalent corners of the BZ, we can write down the eigenvector of the lower band as $`𝐯_+=\frac{1}{\sqrt{2|𝐰|(|𝐰|w_3)}}\left(\begin{array}{c}\hfill w_3|𝐰|\\ \hfill w_1+iw_2\end{array}\right)`$, for $`𝐤V_+`$, and $`𝐯_{}=\frac{1}{\sqrt{2|𝐰|(|𝐰|+w_3)}}\left(\begin{array}{c}\hfill w_1+iw_2\\ \hfill w_3+|𝐰|\end{array}\right)`$, for $`𝐤V_{}`$. There remains some freedom in the division of the BZ. Because $`𝐯_\pm `$ is ill-defined only at $`𝐤_{}^0`$, respectively, we are free to deform this division as long as $`𝐤_{}^0V_\pm `$. This corresponds to the gauge degree of freedom. At $`𝐤V_+V_{}`$, two choices of the wavefunction are different by a phase factor $`e^{i\varphi (𝐤)}𝐯_+=𝐯_{}`$, i.e., $`𝐀_+𝐀_{}=_𝐤\varphi (𝐤)`$, where $`e^{i\varphi (𝐤)}=\frac{w_1iw_2}{\sqrt{w_1^2+w_2^2}}`$. Thus, we get $`\text{Ch}=_{V_+}(𝐀_+𝐀_{})\frac{\text{d}𝐤}{2\pi }=_{V_+}\frac{\text{d}\varphi (𝐤)}{2\pi }=1`$ haldane . The nonzero Chern number implies that the wavefunction cannot be written as a single function for the entire BZ. This also affects the definition of field operators $`a_{(i)𝐤}`$. Second-quantized Hamiltonian is written as $`=_{𝐤,i,j}c_{𝐤j}^{}H_{ji}(𝐤)c_{𝐤i}`$. Field operators $`a_{(i)𝐤}`$ are written as $`a_{(i)𝐤}=_jv_{(i)j}^{}(𝐤)c_{j𝐤}`$. Let $`a_{𝐤\pm }`$ denote the annihilation operators for the lower band when $`𝐤V_\pm `$. Then, $`e^{i\varphi (𝐤)}𝐯_+=𝐯_{}`$ yields $`a_{𝐤+}=e^{i\varphi (𝐤)}a_𝐤`$ at $`𝐤V_+V_{}`$.

We now consider the singlet SC on the lower band. In conventional formalism in the BCS theory of SC, it is assumed that the field operators are continuous in the whole BZ; it is no longer true in the present case. We assume that the Fermi energy crosses only the lower band. The BCS pairing term in the Hamiltonian is

$$\underset{𝐤V_+}{}\mathrm{\Delta }(𝐤)_+a_{𝐤+}^{}a_𝐤^{}+\underset{𝐤V_{}}{}\mathrm{\Delta }(𝐤)_{}a_𝐤^{}a_{𝐤+}^{}.$$
(1)
Singlet SC implies that $`\mathrm{\Delta }(𝐤)_+=\mathrm{\Delta }(𝐤)_{}`$ for $`𝐤V_+`$, i.e. $`\mathrm{\Delta }`$ is an even function of $`𝐤`$. The guiding principle to study topological structure of $`\mathrm{\Delta }(𝐤)`$ is the gauge invariance, corresponding to the freedom in the division of the BZ. It is equivalent to require that the BCS term is continuous at the overlap of two regions; we get $`\mathrm{\Delta }(𝐤)_+=e^{2i\varphi (𝐤)}\mathrm{\Delta }(𝐤)_{}`$ for $`𝐤V_+V_{}`$. The total vorticity $`\nu `$ is

$`\nu ={\displaystyle \underset{m=\pm }{}}{\displaystyle _{V_m}}{\displaystyle \frac{\text{d}\text{arg}\mathrm{\Delta }(𝐤)_m}{2\pi }}={\displaystyle \frac{1}{2\pi }}{\displaystyle _{V_+}}\text{d}\text{arg}\left({\displaystyle \frac{\mathrm{\Delta }_+}{\mathrm{\Delta }_{}}}\right)`$
$`=2{\displaystyle _{V_+}}{\displaystyle \frac{\text{d}\varphi (𝐤)}{2\pi }}=2\text{Ch}=2.`$ (2)
Nonzero $`\nu `$ implies that $`\mathrm{\Delta }(𝐤)`$ can never be expressed as a single continuous function for the entire BZ. To see what the gap function looks like, we take an example of an on-site attraction: $`_U=U_in_in_i`$. We assume that there are two hole pockets surrounding $`𝐤_\pm ^0`$ each. The attractive potential $`V(𝐤,𝐤^{})`$ for the lower band is factorized as a product of a function of $`𝐤`$ and that of $`𝐤^{}`$. We obtain thereby $`\mathrm{\Delta }(𝐤)_\pm =C\frac{w_1(𝐤)iw_2(𝐤)}{|𝐰(𝐤)|}(𝐤V_\pm )`$, where $`C`$ is a complex constant. This has two vortices at $`𝐤=𝐤_\pm ^0`$, with vorticity $`1`$ each. The spin Hall conductivity $`\sigma _{xy}^s`$ spinHall  can be calculated as $`\sigma _{xy}^s=0`$; it is quantized when $`\mathrm{\Delta }(𝐤)0`$ on the FS, and $`\sigma _{xy}^s`$ can be nonzero in general.

We can construct simple 3D models starting from the Haldane’s model. For example, we set $`w_1=t_1_j\mathrm{cos}(𝐤𝐚_j)`$, $`w_2=t_1_j\mathrm{sin}(𝐤𝐚_j)`$, $`w_3=2t_2\mathrm{cos}(k_zc)_j\mathrm{sin}(𝐤𝐛_j)`$ in the Hamiltonian $`H=𝐰𝝈`$, where $`c`$ is the lattice constant in the $`c`$-direction. This is a tight-binding model on a stacked honeycomb lattice with hopping along $`\pm 𝐚_j`$ (nearest-neighbor) and $`\pm 𝐛_j\pm (0,0,c)`$; the latter hopping acquires $`\pi /2`$ phase if it is clockwise in the hexagonal plaquette. This model has $`\rho =_{\alpha ,\beta =\pm }\beta \delta (𝐤𝐤_\alpha ^0+\beta \frac{\pi }{2c}\widehat{𝐜})`$, leading to $`\rho ^\mathrm{\Delta }=2\rho 0`$. In some range of parameters, some of the FS’s encircle a monopole or antimonopole, and the total flux is nonzero for them.

So far we have dealt with noninteracting systems. One may wonder if the topological properties of $`\mathrm{\Delta }(𝐤)`$ discussed above are robust against interactions. In general, robustness of physical phenomena originated from topology is guaranteed by an existence of integer topological numbers, which remain unchanged under an adiabatic change of the Hamiltonian. Although this adiabatic principle usually applies to gapful systems, it applies also to Fermi liquid with interactions. In Fermi liquid, bare electrons turn to quasiparticles which are well-defined near the FS. (Even with interactions, the FS is topologically stable in the momentum space. See Fig. 1 in volovik2 .) The BCS theory is founded on this Fermi liquid theory, and the gap function $`\mathrm{\Delta }(𝐤)`$ is defined in terms of these quasiparticles. The total strength of $`\mathrm{\Delta }(𝐤)`$-monopoles inside the FS is a surface integral of a well-defined function $`\frac{1}{2\pi }_𝐤\times _𝐤\text{arg}\mathrm{\Delta }(𝐤)`$ over the FS, and it remains a well-defined integer even when we include interactions and when $`\mathrm{\Delta }(𝐤)`$ vanishes far from the FS. When we start with noninteracting systems, and adiabatically switch on interactions, this integer remain the same for certain range of the strength of the interactions. We note that also the total strength of monopoles for the wavefunction inside the FS remains an integer in the presence of interactions. Its generalized definition in the presence of interactions is $`\frac{1}{24\pi ^2}ϵ_{\mu \nu \lambda \gamma }\text{Tr}_S𝑑S^\lambda G_{k_\mu }G^1G_{k_\nu }G^1G_{k_\lambda }G^1`$ where $`G(i\omega ,𝐤)`$ is the one-particle Matsubara Green’s function, and $`S`$ is the three-dimensional surface surrounding the FS in the four-dimensional frequency-momentum space $`(\omega ,𝐤)`$ (see Eq. (65) in volovik2 ). It is a well-defined integer, and topologically stable volovik2 . After applying the Stokes’ theorem, the integrand is finite only on the FS, and the integral is equal to the number of monopoles inside the FS. In 3D, Fermi liquid picture is valid even with interactions, and modification of Fermi velocity $`v_F`$ or of quasiparticle residue $`Z`$ does not change the number of monopoles inside the FS. Thus the total monopole strength remains unchanged. It changes only when the topology of the FS’s changes.

Let us discuss on applicability of our theory to real materials. As mentioned above, we should search through magnetic SC with (a) equal-spin pairing without inversion symmetry or (b) opposite-spin pairing. The case (a) is likely in ferromagnetic SC, while (b) is likely in antiferromagnetic SC. Band structure calculations have been concentrating on energy eigenvalues on high-symmetry directions. Almost no attention has been paid to existence of band-touching points, or to phases of wavefunctions. Thus, we do not have at present sufficient information to see which materials are candidates for our theory. In order to consider how often nontrivial topological structure appears, we should instead resort to simplified models. In onoda , the authors consider a ferromagnetic tight-binding model on the 2D square lattice with $`t_{2g}`$ orbitals, with physically reasonable values of parameters. The gauge flux $`B(𝐤)_z`$ has sharp peaks, where the bands approach each other in energy. If we go to the 3D by adding an extra dimension along $`k_z`$, the bands will touch each other at some points, because band-touching points (diabolical points) has codimension three volovik ; volovik2 . Actually in shindou , the diabolical points, i.e. the (anti)monopoles, are found in a model for an antiferromagnet. These points are topologically stable. Thus nontrivial topology of the wavefunctions and the gap functions, as discussed in the present Letter, is expected in many materials.

We can argue exciton condensation in the similar way. In 2D, total vorticity of the order parameter is equal to the difference of Chern numbers of the bands involved. In 3D, the monopole density of the order parameter is the difference of those for wavefunctions of the two bands.

In conclusion, we considered the SC in a band with nontrivial topology. It is formulated in terms of a gauge flux distribution in the $`𝐤`$-space. Although the gap function $`\mathrm{\Delta }(𝐤)`$ depends on a detail of attractive potential, topological structure of $`\mathrm{\Delta }(𝐤)`$ is determined solely from that of the normal-state wavefunction. The effect of disorder is an interesting issue, and is left for future studies.

###### Acknowledgements.

The authors thank helpful discussion with R. Shindou and M. Sigrist. We acknowledge support by Grant-in-Aids from the Ministry of Education, Culture, Sports, Science and Technology of Japan.


text:
# 1 Introduction

## 1 Introduction

When a two-level atom and a photon mode undergo a Jaynes–Cummings interaction , they become entangled as a rule. Various aspects of atom-field entanglement of this kind have been studied in past years, whereby the emphasis is usually — but not always  — on the consequent atom-atom entanglement that results from the interaction of several atoms with the same photon mode.

Clearly, if two atoms interact with the photons at the same time, rather strong atom-atom entanglement is expected, and this is indeed the case (see, e.g., the recent paper by Masiak ). For similar reasons one also gets strong entanglement if particular initial internal states of the atoms are suitably combined with well chosen interaction times (see, e.g., ).

Much less obvious are the properties of the atom-atom entanglement that is generated in standard one-atom maser (OAM) experiments . In general terms it is known that there are correlations among emerging atoms, that they manifest themselves in non-Poissonian counting statistics , for instance, and that measurements on consecutive atoms can be inconsistent with a Bell-type inequality .

But a systematic study of the properties of the joint state of two atoms that have emerged from the resonator of a OAM experiment has not been performed beyond the extraction of some correlation data of experimental relevance . It is the objective of the present contribution, which is a progress report in character, to contribute to this systematic study.

We shall here focus on what is perhaps the simplest situation in this context, namely two atoms where the second traverses the resonator very shortly after the first leaves. One anticipates correctly that such atom pairs exhibit particularly strong entanglement because the dissipative degradation of the photon field can be disregarded for the brief period between the two atoms.

## 2 Successive atoms in a one-atom maser

To set the stage, we first recall some basics of OAM theory (see, e.g., ). The resonant Jaynes–Cummings dynamics evolves the joint state of a two-level atom (excited state: $`|𝖾`$, ground state: $`|𝗀`$) and the photon field (state of definite photon number: $`|n`$, with $`n=0,1,2,\mathrm{}`$) in accordance with

$`|𝖾,n,\mathrm{}`$ $``$ $`𝐂|𝖾,n,\mathrm{}i𝐒|𝗀,n,\mathrm{}`$ (1)
$`=`$ $`|𝖾,n,\mathrm{}\mathrm{cos}\left(\phi \sqrt{n+1}\right)i|𝗀,n+1,\mathrm{}\mathrm{sin}\left(\phi \sqrt{n+1}\right),`$
where the ellipsis indicates quantum numbers that are not affected (such as those referring to other atoms) and

$$𝐂=\mathrm{cos}\left(\phi \sqrt{aa^{}}\right),𝐒=a^{}\frac{\mathrm{sin}\left(\phi \sqrt{aa^{}}\right)}{\sqrt{aa^{}}}$$
(2)
are essentially trigonometric functions of the photon ladder operators $`a`$ and $`a^{}`$. The accumulated Rabi angle $`\phi `$ is the net measure of the interaction strength and duration.

Prior to the arrival of the first atom, the photon field is in the steady state, specified by the statistical operator $`\varrho ^{(\mathrm{ss})}`$ that is given by

$$\varrho ^{(\mathrm{ss})}(a^{}a)=\varrho ^{(\mathrm{ss})}(0)\underset{n=1}{\overset{a^{}a}{}}\left[\frac{\nu }{\nu +1}+\frac{N_{\mathrm{ex}}}{\nu +1}\frac{\mathrm{sin}^2\left(\phi \sqrt{n}\right)}{n}\right].$$
(3)
Here, the normalized pump rate $`N_{\mathrm{ex}}`$ is the mean number of atoms traversing the resonator in one photon lifetime, $`\nu `$ is the number of thermal photons in the cavity field, and the value of $`\varrho ^{(\mathrm{ss})}(0)`$ is determined by the normalization of $`\varrho ^{(\mathrm{ss})}`$ to unit trace.

The final state of the first atom is thus given by

$$\varrho ^{(1\mathrm{s}\mathrm{t})}=|𝖾\mathrm{tr}\left\{𝐂\varrho ^{(\mathrm{ss})}𝐂\right\}𝖾|+|𝗀\mathrm{tr}\left\{𝐒\varrho ^{(\mathrm{ss})}𝐒^{}\right\}𝗀|=\frac{1}{2}(1+\stackrel{}{s}\stackrel{}{\sigma }),$$
(4)
where

$$\stackrel{}{s}=(0,0,s)\text{with}s=\mathrm{tr}\left\{(𝐂^2𝐒^{}𝐒)\varrho ^{(\mathrm{ss})}\right\}$$
(5)
is the Bloch vector for the first atom and

$$\stackrel{}{\sigma }=(\sigma _x,\sigma _y,\sigma _z)=(|𝗀𝖾|+|𝖾𝗀|,i|𝗀𝖾\left|i\right|𝖾𝗀|,|𝖾𝖾||𝗀𝗀|)$$
(6)
is the Pauli vector operator of the familiar pseudo-spin formalism for two-level atoms.

For the second atom, we denote the Pauli vector operator by $`\stackrel{}{\tau }`$ and the numerical Bloch vector by $`\stackrel{}{t}`$. Since the second atom is assumed to arrive immediately after the first has left — this is the “perhaps simplest situation” alluded to in the Introduction — its statistical operator is given by

$$\varrho ^{(2\mathrm{n}\mathrm{d})}=\frac{1}{2}(1+\stackrel{}{t}\stackrel{}{\tau })$$
(7)
with

$$\stackrel{}{t}=(0,0,t),t=\mathrm{tr}\left\{(𝐂^2𝐒^{}𝐒)(𝐂\varrho ^{(\mathrm{ss})}𝐂+𝐒\varrho ^{(\mathrm{ss})}𝐒^{})\right\}.$$
(8)
Rare exceptions aside, the values of $`s`$ and $`t`$ will be different, which is, of course, a sign that the cavity field “remembers” that the first atom has just passed through the resonator. Photon dissipation makes this memory fade in the course of time and, therefore, it is quite obvious that the correlations between two successive atoms are strongest if the second comes immediately after the first. In the present paper, we are dealing with this particularly favorable situation. Results about more general cases will be reported elsewhere.

Now turning to the correlated statistics of two successive atoms, we note that, in addition to the Bloch vectors $`\stackrel{}{s}`$ and $`\stackrel{}{t}`$, the statistical operator $`\varrho ^{(\mathrm{both})}`$ of the joint two-atom state,

$$\varrho ^{(\mathrm{both})}=\frac{1}{4}\left(1+\stackrel{}{s}\stackrel{}{\sigma }+\stackrel{}{t}\stackrel{}{\tau }+\stackrel{}{\sigma }\genfrac{}{}{0pt}{}{}{C}\stackrel{}{\tau }\right),$$
(9)
involves the cross dyadic $`\genfrac{}{}{0pt}{}{}{C}`$, which happens to be diagonal,

$$\genfrac{}{}{0pt}{}{}{C}=\left(\begin{array}{ccc}u& 0& 0\\ 0& u& 0\\ 0& 0& v\end{array}\right)$$
(10)
with

$`u`$ $`=`$ $`\mathrm{tr}\left\{(𝐒^{}\mathrm{𝐂𝐒𝐂}+\mathrm{𝐂𝐒}^{}\mathrm{𝐂𝐒})\varrho ^{(\mathrm{ss})}\right\}`$
$`\text{and}v`$ $`=`$ $`\mathrm{tr}\left\{(𝐂^2𝐒^{}𝐒)(𝐂\varrho ^{(\mathrm{ss})}𝐂𝐒\varrho ^{(\mathrm{ss})}𝐒^{})\right\}.`$ (11)
By construction, all eigenvalues of $`\varrho ^{(\mathrm{both})}`$ are non-negative, which is to say that the four parameters $`s`$, $`t`$, $`u`$, and $`v`$ obey the inequalities

$$1v\sqrt{4u^2+(st)^2},1+v|s+t|.$$
(12)
The difference between $`\varrho ^{(\mathrm{both})}`$ and the two-atom product state $`\varrho ^{(1\mathrm{s}\mathrm{t})}\varrho ^{(2\mathrm{n}\mathrm{d})}`$,

$$\mathrm{\Delta }\varrho =\varrho ^{(\mathrm{both})}\varrho ^{(1\mathrm{s}\mathrm{t})}\varrho ^{(2\mathrm{n}\mathrm{d})}=\frac{1}{4}\stackrel{}{\sigma }\genfrac{}{}{0pt}{}{}{E}\stackrel{}{\tau },$$
(13)
is measured in terms of the entanglement dyadic $`\genfrac{}{}{0pt}{}{}{E}=\genfrac{}{}{0pt}{}{}{C}\stackrel{}{s}\stackrel{}{t}`$. Here its only non-zero elements are the values $`u,u,vst`$ on the diagonal. Their absolute values are the characteristic values of $`\genfrac{}{}{0pt}{}{}{E}`$, and the largest of them is the maximal “degree of correlation” , which is one measure of the strength of the entanglement between the two atoms.

Another and perhaps more telling measure is the absolute size of $`\mathrm{\Delta }\varrho `$, its trace norm, $`\mathrm{tr}\left\{\left|\mathrm{\Delta }\varrho \right|\right\}`$, the sum of the moduli of the eigenvalues of $`\mathrm{\Delta }\varrho `$. These four eigenvalues are $`\frac{1}{4}(2uvst),\frac{1}{4}(2uvst)`$, and twice $`\frac{1}{4}(vst)`$, so that we get

$$\mathrm{tr}\left\{\left|\mathrm{\Delta }\varrho \right|\right\}=\{\begin{array}{cc}|vst|\hfill & \text{if }2|u||vst|\text{,}\hfill \\ \frac{1}{2}|vst|+|u|\hfill & \text{if }2|u||vst|\text{,}\hfill \end{array}$$
(14)
for this measure of entanglement.

## 3 Degree of separability

Although an entangled two-atom state $`\varrho ^{(\mathrm{both})}`$ is, by definition, not equal to the product $`\varrho ^{(1\mathrm{s}\mathrm{t})}\varrho ^{(2\mathrm{n}\mathrm{d})}`$ of the two single-atom states contained in it, it may very well be a convex sum of such products. Then one speaks of a *separable* state,

$$\varrho ^{(\mathrm{both})}\text{ is separable if}\varrho ^{(\mathrm{both})}=\underset{k}{}w_k\varrho _k^{(1\mathrm{s}\mathrm{t})}\varrho _k^{(2\mathrm{n}\mathrm{d})}\text{with }w_k>0.$$
(15)
Entangled, but separable states exhibit correlations that can be modeled classically: It is as if we had one of the products $`\varrho _k^{(1\mathrm{s}\mathrm{t})}\varrho _k^{(2\mathrm{n}\mathrm{d})}`$ without knowing which one, the $`k`$-th occurring with probability $`w_k`$.

The Peres–Horodeccy criterion  makes it easy to check whether the $`\varrho ^{(\mathrm{both})}`$ of (9)–(14) is separable for the given set of parameter values. One finds that

$$\varrho ^{(\mathrm{both})}\text{ is separable if}1+v\sqrt{4u^2+(s+t)^2}\text{and only then}.$$
(16)
In particular, we have a separable $`\varrho ^{(\mathrm{both})}`$ whenever $`u=0`$.

According to what is said in the context of (15), truly non-classical correlations require a non-separable state. Now, there are non-separable states that differ very little from separable ones and others that differ by much. A quantitative measure for a non-separable state’s proximity to separable ones is the *degree of separability* $`𝒮`$, which is defined by the Lewenstein–Sanpera decomposition of $`\varrho ^{(\mathrm{both})}`$ ,

$$\varrho ^{(\mathrm{both})}=𝒮\varrho _{\mathrm{sep}}^{(\mathrm{both})}+(1𝒮)\varrho _{\mathrm{pure}}^{(\mathrm{both})},$$
(17)
where $`\varrho _{\mathrm{sep}}^{(\mathrm{both})}`$ is separable, $`\varrho _{\mathrm{pure}}^{(\mathrm{both})}`$ is a pure state, and $`𝒮`$ is maximal. The smaller the value of $`𝒮`$, the stronger are the non-classical correlations in $`\varrho ^{(\mathrm{both})}`$.

For states of the particularly simple 4-parametric structure (9)–(14), one can use known methods  to find the decomposition (17), with this outcome: The pure part is of the form

$$\varrho _{\mathrm{pure}}^{(\mathrm{both})}=\frac{1}{4}\left[1+p\sigma _zp\tau _z+q\left(\sigma _x\tau _x+\sigma _y\tau _y\right)\sigma _z\tau _z\right]$$
(18)
with $`q=\mathrm{sgn}(u)\sqrt{1p^2}`$ and the value of $`p`$ is determined as follows. We have

$$p=0\text{if}\mathrm{\Lambda }^2(1+v)\mathrm{\Lambda }+st0$$
(19)
with

$$\mathrm{\Lambda }=1|u|+\frac{1}{2}\sqrt{(1+v)^2(s+t)^2}$$
(20)
and get $`p`$ as the solution of

$$\frac{1\mathrm{\Lambda }}{\sqrt{1p^2}}\left[(1v)(st)p\right]=(1\mathrm{\Lambda })^2v+st$$
(21)
when (19) does not apply. The degree of separability is then given by

$$𝒮=1\frac{1\mathrm{\Lambda }}{\sqrt{1p^2}},$$
(22)
so that $`𝒮=\mathrm{\Lambda }`$ when $`p=0`$, and $`𝒮<\mathrm{\Lambda }`$ when $`p0`$.

## 4 Examples and summary

In Fig. 1 we show some numerical results obtained for the OAM parameters $`(N_{\mathrm{ex}},\nu )=(1,0)`$, $`(1,0.2)`$, $`(3,0)`$, and $`(5,0)`$. All are for rather small fluxes, namely $`N_{\mathrm{ex}}=1`$ in plots (a) and (b), $`N_{\mathrm{ex}}=3`$ in plot (c), and $`N_{\mathrm{ex}}=5`$ in plot (d). In the present context, small fluxes have the advantage that the steady state $`\varrho ^{(\mathrm{ss})}`$ does not have a large average photon number, so that the additional photon that can be deposited by the first atom may lead to a substantial change in the cavity field. The field then encountered by the second atom is markedly different from $`\varrho ^{(\mathrm{ss})}`$ and, therefore, there is a good chance to generate a strong entanglement between the atoms.

Put differently, with a large flux, $`N_{\mathrm{ex}}1`$, one will have a steady state with a photon number distribution spread over a wide range. The statistical average over the various values of $`n`$ in (1) will then tend to degrade the atom-atom correlations. Indeed, the comparison of plots (a), (c), and (d) — where $`N_{\mathrm{ex}}=1`$, $`3`$, and $`5`$, respectively — demonstrates that both $`\mathrm{tr}\left\{|\mathrm{\Delta }\varrho |\right\}`$ and $`1𝒮`$ have somewhat smaller values for larger $`N_{\mathrm{ex}}`$. The effect on $`𝒮`$ is particularly pronounced as the total dark area is much smaller in plot (c) than in plot (a), and even smaller in plot (d).

Similarly, one expects that the noise associated with thermal photons has a negative influence on the atom-atom correlations. This is confirmed by the comparison between plots (a) and (b), where $`\nu =0`$ and $`\nu =0.2`$, respectively. Again we notice some reduction in $`\mathrm{tr}\left\{|\mathrm{\Delta }\varrho |\right\}`$ and a marked loss of dark area.

In summary, then, we have demonstrated that a one-atom maser in standard operation can give rise to rather strong entanglement between successive atoms, provided the parameters are chosen judiciously (low flux, low temperature). For some parameter ranges, one obtains non-separable two-atom states with truly non-classical correlations.

We have confined the discussion to the simple situation of atoms that traverse the resonator in immediate succession. Clearly, one must allow for a more reasonable time gap when aiming at a comparison with real-life experimental data. Work in this direction is in progress and results will be reported in due course. Further, the entanglement properties of chains of more than two atoms are not understood to a satisfactory extent. So far there have only been studies of bulk quantities, such as the mean number of successive atoms in the same final state , and surely there is a lot of finer detail that deserves more attention.

## Acknowledgments

BGE’s work is supported in part by a TITF initiative of Texas A&M University.

## Figures


text:
# References

Optical soliton solutions of the quintic complex Swift-Hohenberg equation
A. Ankiewicz1, K. Maruno2 and N. N. Akhmediev3

1 Applied Photonics Group, and 3 Optical Sciences Centre, Research School of Physical Sciences and Engineering,

The Australian National University, Canberra ACT 0200, Australia

2 Research Institute for Applied Mechanics, Kyushu University,Kasuga, Fukuoka, Japan.

1. Introduction. Complicated optical systems with gain and loss can be described by the complex Swift-Hohenberg (S-H) equation. These include synchronously pumped optical parametric oscillators, pattern formation in convection cells and passively mode-locked lasers with special spectral features -. The main difference between the complex S-H equation and previously-studied master equations lies in its more involved spectral filtering term. At the same time, these complications do not allow us to analyse the solutions easily. In fact, it was not clear that such solutions may exist at all -. In this work, we study the quintic complex S-H equation in 1D and report various new exact solutions.

The normalized quintic complex S-H equation is:

$$i\psi _z+\frac{D}{2}\psi _{tt}+|\psi |^2\psi +(h+is)\psi _{tttt}+(\nu i\mu )|\psi |^4\psi =i\delta \psi +i\beta \psi _{tt}+iϵ|\psi |^2\psi .$$

In mode-locked laser applications, $`z`$ is the propagation distance or the cavity round-trip number (treated as a continuous variable), $`t`$ is the retarded time, $`D`$ is the 2nd order dispersion, $`h`$ is the 4th order dispersion, $`ϵ`$ is a nonlinear gain (or 2-photon absorption if negative) and $`\delta `$ (usu. negative) represents a constant gain or loss. The band-limited gain (e.g. due to an EDFA, where the gain band may be about 30 nm around 1.5 microns) is represented by $`\beta `$ (parabolic spectrum shape) and $`s`$ (4th order correction). We find exact forms for a range of solitons, including bright and dark cases, and both chirped and unchirped forms. They have various features which differentiate them from solitons of the complex Ginzburg-Landau equation . These first solutions can give some clues for analysing more involved solutions.

2. Analysis. Let $`\psi =f(t)exp[i\mathrm{\Omega }z]`$. This reduces the quintic complex S-H equation to an ordinary differential eqn. in $`f`$:

$$\mathrm{\Omega }f+(\frac{D}{2}i\beta )f^{\prime \prime }(t)+(1iϵ)|f|^2fi\delta f+(h+is)f^{\prime \prime \prime \prime }(t)+(\nu i\mu )|f|^4f=0.$$
(1)
We will show that basic \[unchirped\] bright and dark solitons exist, and that there are also chirped bright and dark solitons.

3.1 Bright soliton If we take $`f=cgsech(gt)`$, then we note that $`f^{\prime \prime }(t)/[g^2f(t)]=\mathrm{\hspace{0.17em}\hspace{0.17em}2}tanh^2(gt)1`$ and that $`f^{\prime \prime \prime \prime }(t)/[g^4f(t)]=\mathrm{\hspace{0.17em}\hspace{0.17em}5}28tanh^2(gt)+24tanh^4(gt).`$ On dividing through by $`f`$, it is clear that the coefficients of $`tanh^n(gz)`$ for $`n=0,2,4`$ can be set to zero to reduce the problem to an algebraic one and obtain the solution. From the $`n=4`$ term , we find a consistency condition on the equation parameters: $`\nu =h\mu /s.`$ Thus $`c^4=\mathrm{\hspace{0.17em}24}s/\mu (=24h/\nu )`$. From the other terms we find $`g,c,etc`$ and conditions we need to impose as constraints on the equation parameters. Thence we can easily find $`\beta ,\delta ,g`$ and $`\mathrm{\Omega }`$ Note that , if $`ϵ=0`$ the solution simplifies, as then $`g=\frac{\beta }{10s}`$, so

$$f=\sqrt{\frac{\beta }{5}}(\frac{6}{s\mu })^{\frac{1}{4}}sech[\sqrt{\frac{\beta }{10s}}t].$$

3.2 Black soliton. Here we take If we take $`f=cgtanh(gt)`$, then we see that $`f^{\prime \prime }(t)/[g^2f(t)]=2[tanh^2(gt)1],`$ and that $`f^{\prime \prime \prime \prime }(t)/[g^4f(t)]=8[tanh^2(gt)1][3tanh^2(gt)2].`$ Again, on dividing by $`f`$, it is clear that the coefficients of $`tanh^n(gt)`$ for $`n=0,2,4`$ can be set to zero to reduce the problem, and we then determine all the solution parameters .

4.1. Chirped bright soliton. This involves the function $`f=a(t)exp[idloga(t)]`$, where $`a(t)=gcsech(gt)`$. With this form, the derivatives can still be written in a convenient way. For example, $`f^{\prime \prime }(t)/[g^2f(t)]=(di)[i+(d2i)tanh^2(gt)].`$

We need to find the roots of a 4th order polynomial in $`d`$:

$$(\mu h+s\nu )(d^435d^2+24)+\mathrm{\hspace{0.17em}10}d(\nu hs\mu )(5d^2)=\mathrm{\hspace{0.17em}0}.$$
(2)
Thus $`d`$ depends on a balance of the highest order derivative (4th) and strongest nonlinearity and is quite different from that of the CGLE  or CGLE with an integral term  where $`d`$ is determined by $`\beta `$ and $`ϵ`$. Eqn(2) provides insight into the chirpless ($`d=0`$) case, as we see that $`d=0`$ is a root of eqn.(2) when $`\mu h+s\nu =0`$.

4.2 Chirped black soliton Here $`f(t)=a(t)exp[idlogb(t)],`$ where $`a(t)=rgtanh(gt)`$ and $`b(t)=gsech(gt)`$, and we proceed as above to find the solution.

5. Energy and momentum balance. As any other equation describing dissipative systems, the quintic complex S-H equation does not have any conserved quantities. Instead, we can write balance equations for the energy and momentum. A study of the quintic complex S-H eqn., using the energy balance approach (,) leads to the following evolution equation for the energy:

$$\frac{d}{dz}\underset{\mathrm{}}{\overset{\mathrm{}}{}}|\psi |^2𝑑t=2\underset{\mathrm{}}{\overset{\mathrm{}}{}}\left[\delta |\psi |^2\beta \left|\frac{\psi }{t}\right|^2s|\psi _{tt}|^2+ϵ|\psi |^4+\mu |\psi |^6\right]𝑑t.$$
(3)
By definition, the right hand side of this equation is the rate of change of the energy. The term with $`s`$ here is new in comparison with similar equation for CGLE . The form of this term reflects its role as a higher-order band-limited gain. For any stationary exact solution, we need the r.h.s. of eqn.(3), and also the rate of change of the momentum, to be zero. This provides a way of finding or checking solutions. The balance equation is also an important tool in studying the interaction between the pulses .

Conclusion. Solutions presented here are novel examples of exact solutions which exist for the quintic complex S-H equation. As for the Ginzburg -Landau equation  they certainly do not cover the whole set of possible solutions. In fact, they only represent a small subset of a variety of soliton-like solutions. Other solutions have to be studied numerically and this may take considerable amount of simulations. However, finding the exact solutions is an important step in analysing the laser system with more involved spectral properties which is described by the quintic complex S-H equation.