Shotaro Kikuchihara

Masami Ouchi

Yoshiaki Ono

Ken Mawatari

Jacopo Chevallard

Yuichi Harikane

Takashi Kojima

Masamune Oguri

Gustavo Bruzual

Stephane Charlot

Early Low-Mass Galaxies and Star-Cluster

Keywords: galaxies, formation -galaxies, high-redshift -galaxies, luminosity function, mass function -globular clusters, general -gravitational lensing, strong

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

skiku@icrr.u-tokyo.ac.jp

The stellar mass of a galaxy is an indicator of the masses that are aggregated through the previous star formation and merging processes. Various observations and simulations have found that the stellar masses correlate with basic properties of galaxies, such as the star formation rate (SFR; e.g., McLure et al. 2011;Speagle et al. 2014;Salmon et al. 2015;Santini et al. 2017;Iyer et al. 2018), age (e.g., Sparre et al. 2015), metallicity (e.g., Kojima et al. 2017;Barber et al. 2018), and size (Trujillo et al. 2004;van der Wel et al. 2014;Lange et al. 2015). These facts imply that the stellar mass plays a critical role in understanding the formation and evolution of galaxies in the early universe, which still are major open questions in astronomy today.

The number density of galaxies per stellar-mass (M ) interval, i.e., the galaxy stellar mass function (GSMF), is often used to study the evolution of the total stellar mass in the universe. The low-mass ends of the GSMFs are especially worth investigating at high redshift, because low-mass galaxies are expected to be dominant in the early universe according to the hierarchical cosmology. It is also worth mentioning that the low-mass galaxies (M ∼ 10 6 -10 7 M ) at z ∼ 6 -9 are predicted to grow up to have stellar masses comparable to that of the Milky Way at z ∼ 0 (e.g., Behroozi et al. 2013Behroozi et al. , 2018)).

However, the low-mass ends of the GSMFs at high redshift are poorly constrained due to the lack of the sample of high-redshift low-mass galaxies caused by insufficient depths of observations. These observational limits can be resolved by the capability of the Wide Field Camera 3 (WFC3/IR) on board the Hubble Space Telescope (HST) together with the strong gravitational lensing effect caused by galaxy clusters. In this manner, the Hubble Frontier Fields (HFF; Coe et al. 2015;Lotz et al. 2017) spent more than 800 orbits of the HST to survey six galaxy cluster fields, extending the faint end of the rest-frame ultraviolet (UV) luminosity functions (LFs) to a intrinsic UV magnitude M UV ≈ -14 mag at z 6 (Finkelstein et al. 2015;Laporte et al. 2016;McLeod et al. 2016;Livermore et al. 2017;Ishigaki et al. 2018;Atek et al. 2018;Bhatawdekar et al. 2018). Here, one should study the GSMFs complementary to the UV LFs that represent the fundamental physical quantity of the stellar mass.

The lensing effect has brought further benefit to studies of galaxy structures and morphologies. Lensed galaxies are stretched along critical curves, allowing the structures to be studied at high spatial resolution. Note that this can be fully exploited with the high resolving power of the WFC3/IR. Owing to magnification by the lensing effects, understanding the size evolution has been significantly advanced by recent studies (e.g., Shibuya et al. 2015;Bouwens et al. 2017;[START_REF] Kawamata | [END_REF]. Bouwens et al. (2017) compare the effective radii R e and the stellar masses M of the z ∼ 6 -8 galaxies to study the analogy of the high-redshift galaxies to the local stellar systems, reporting that some galaxies have the values of R e and M comparable to those of the globular clusters (GCs), super star complexes, or star cluster complexes in the local universe. Here, Bouwens et al. (2017) assume M -M UV relations, i.e. stellar mass to UV luminosity L UV ratios M /L UV , to estimate M with M UV measurements. In Bouwens et al. (2017), the assumed M -M UV relations are based on two representative models with star-formation duration time of 10 and 100 Myr (including the time difference by a factor of ∼ 4) that have not been tested with observational data yet. The M -M UV relation should be determined by observations with no such assumptions, and the R e -M relation should be studied.

The structure of this paper is as follows. We describe the data and sample in detail in Section 2. The M -M UV relations are obtained by the staking analysis and the spectral energy distribution (SED) fitting technique in Section 3. In Section 4, we show the results and the discussions on the GSMFs and the R e -M relations of the galaxies. Finally, we give a summary in Section 5.

Throughout this paper, we adopt a cosmology with Ω m,0 = 0.3, Ω Λ,0 = 0.7, and H 0 = 70 km s -1 Mpc -3 . Magnitudes are given in the AB system (Oke & Gunn 1983). Densities and sizes of the galaxies are measured in comoving and physical scales, respectively. We adopt the Chabrier (2003) initial mass function (IMF) in a mass range of 0.1 -100 M to estimate stellar masses. In this paper, all of the stellar masses taken from the previous studies are converted to those estimated with the Chabrier (2003) IMF.

DATA AND SAMPLES

HST Data and Samples

In our analysis, we make use of the image mosaics obtained in the HFF program, which targets six cluster fields-Abell 2744, MACS J0416.1-2403, MACS J0717.5+3745, MACS J1149.6+2223, Abell S1063, and Abell 370 (hereafter A2744, M0416, M0717, M1149, A1063, and A370, respectively)-and their accompanying six parallel fields. All of the 12 fields were observed with the three bands of the Advanced Camera for Surveys (ACS) and four bands of the WFC3/IR; F435W (B 435 ), F606W (V 606 ), F814W (i 814 ), F105W (Y 105 ), F125W (J 125 ), F140W (JH 140 ), and F160W (H 160 ). We utilize the drizzled and weight images that were produced by Shipley et al. (2018) in the manner summarized below. 1 First, the HFF v1.0 images were downloaded from the MAST archive. 2 The point-spread functions (PSFs) of these images were homogenized to those of the H 160 images. The PSF FWHM of the homogenized images is ≈ 0. 18. Second, the bright cluster galaxies (bCGs) were modeled and subtracted from the images to avoid the diffuse intracluster light (ICL) in photometry of the background faint sources. All of the images have a pixel scale of 0. 06. We correct for the Galactic extinction using the values given by Shipley et al. (2018).

We divide each image into 3×3 grid cells, and measure the limiting magnitude in each cell (∼ 1 arcmin 2 ). This is because limiting magnitudes are not homogeneous due to the ICL (e.g., Montes & Trujillo 2014;Ishigaki et al. 2015;Kawamata et al. 2016). The 5σ limiting magnitudes in the H 160 band images are ≈ 28.4 -29.2 mag in a 0. 35-diameter circular aperture.

We use the galaxy sample selected by [START_REF] Kawamata | [END_REF], which consists of 350 i-, 64 Y -, and 39 YJdropouts (or z ∼ 6 -7, z ∼ 8, z ∼ 9 Lyman break galaxies, respectively). Photometry is reperformed in the following way. First, we measure the aperture magnitude m ap with a diameter of D ap ≡ 0. 35 at the position of the dropouts, using the IRAF task phot [START_REF] Tody | Proc. SPIE[END_REF][START_REF] Tody | Astronomical Data Analysis Software and Systems II[END_REF]. Second, we apply an aperture correction. To evaluate the aperture correction term, we create a median-stacked J 125 -band image of the i-dropouts with PSF homogenization, and measure the aperture flux of the stacked dropout with changing the aperture diameter. We define D 99 by the diameter of the aperture that includes ≥ 99 % of the total flux is included. The aperture correction term, c ap , is defined as c ap ≡ m ap (D ap )m ap (D 99 ), where m ap (d) is a aperture magnitude for a given diameter d. For B 435 and V 606 bands, c ap is defined as the same value defined in the i 814 band, because c ap cannot be estimated due to low signal-to-noise (S/N) ratios in the B 435 and V 606 bands. The values of c ap are 0.98, 0.98, 0.98, 0.93, 0.95, 0.94, and 0.90 magnitudes in the B 435 , V 606 , i 814 , Y 105 , J 125 , JH 140 , and H 160 bands, respectively. Finally, we estimate the total magnitude m tot as m tot = m ap (D ap )c ap .

In our photometry, some dropouts may be detected at the bands blueward of the Lyman break. We remove the dropouts with a ≥ 2σ detection either in B 435 or V 606 band for i-dropouts, and in B 435 , V 606 , or i 814 band for Y -and YJ-dropouts. For our dropouts, we apply criteria of ≥ 3σ detections in J 125 , JH 140 , and H 160 bands for i-, Y -, and YJ-dropouts, respectively. Our 2 http://www.stsci.edu/hst/campaigns/frontier-fields/ final sample consists of 267 i-dropouts, 54 Y -dropouts, and 36 YJ-dropouts (357 in total).

VLT and Keck Data

We use the K s band images obtained by the Kband Imaging of the Frontier Fields (KIFF) program (Brammer et al. 2016). The deep K s images of the VLT/HAWK-I (Keck-I/MOSFIRE) are available in both the cluster and the parallel fields of A2744, M0416, A1063, and A370 (M0717 and M1149). The HAWK-I and MOSFIRE images have the PSF FWHM of ≈ 0. 4 -0. 5. We utilize the drizzled and weight images of Shipley et al. (2018), whose details are described in Section 2.1. The pixel scale of the images is 0. 06 that matches to that of the HST images. Photometry and limiting-magnitude measurements are conducted in the same manner as those in Section 2.1, but with D ap = 0. 6 and c ap = 0.98. The 5σ limiting magnitudes are ≈ 25.3 -26.2 mag.

Spitzer Data

We utilize the drizzled and weight images of ch1 (3.6 µm) and ch2 (4.5 µm) of the IRAC. The photometric data that were taken by 2016 December are combined by Shipley et al. (2018), who also applies bCG subtraction and Galactic extinction correction. In both ch1 and ch2, the PSF FWHM is ≈ 1. 7 -2. 0, and pixel scales are 0. 3. Photometry and limiting-magnitude measurements are conducted in the same manner as described in Section 2.1, but with D ap = 3. 0 for both ch1 and ch2, and c ap = 0.52, 0.55 for ch1 and ch2, respectively (Ono et al. 2010). The 5σ limiting magnitudes are ≈ 24.5 -25.5 mag.

Lens Models

We apply the best-fit magnification factors µ to estimate the intrinsic magnitudes and stellar masses of the dropouts. The magnification factors that we use are calculated with the glafic (Oguri 2010) parametric models by Kawamata et al. (2016[START_REF] Kawamata | [END_REF]. 3 Priewe et al. (2017) estimate the uncertainties of the eight mass models of M0416, using a fractional normalized median absolute deviation, which is defined as

fNMAD ≡ 1.4826 × median|µ m (θ) -μ|/μ.
(1)

Here µ m (θ) is the magnification factor calculated with a given model m at a given position θ, and μ is the median of the magnification factors. It is found that fNMAD ≈ 0.3 (0.7) at μ ≈ 2 (40) in M0416, which implies that magnification factors of different models differ only by a factor of ≈ 2. The magnification uncertainties do not change our conclusions (see also Ishigaki et al. 2015;Meneghetti et al. 2017;[START_REF] Kawamata | [END_REF].

STELLAR POPULATION OF THE HIGH-Z GALAXIES

Stacking Analysis

To create images with high S/N ratios, we conduct a stacking analysis in the following way. We take J 125 , JH 140 , and H 160 band magnitudes as the rest-frame UV apparent magnitudes m UV of the i-, Y -, and YJdropouts, respectively. The intrinsic absolute magnitude of a dropout is calculated via

M UV = m UV -5 log d L (z) 10 pc +2.5 log[µ(z)]+2.5 log(1+z),
(2) where d L (z) is the luminosity distance to the dropout, and µ(z) is the best-fit magnification factor at the position of the dropout (see Section 2.4). We apply z = 6, 8, and 9 for the i-, Y -, and YJ-dropouts, respectively, for consistency with the assumption in the size measurements [START_REF] Kawamata | [END_REF]. We divide our dropout samples into subsamples by the values of M UV at z ∼ 6 -7, 8, and 9. A summary of the subsamples is shown in Table 1.

In each band, we cut out images centered at the positions of dropouts, and divide the pixel counts by µ of the dropouts. Here we assume that the values of µ do not change within sizes of the dropout galaxies, because the sizes of the dropout galaxies are too small to significantly impact on the dropout galaxy fluxes. We then median-stack the images of the dropouts for each subsample with iraf task imcombine.

Figures 1,2, and 3 show the stacked images of the z ∼ 6 -7, 8, and 9 subsamples, respectively. Magnitudes for the subsamples are measured in the same manner as those for the individual dropouts in Section 2.1. To estimate uncertainties of the total fluxes (or the total magnitudes), we first follow the three steps; 1) randomly selecting positions in the sky area of the grid cell where the subsamples are located, 2) generating medianstacked sky noise images, and 3) performing aperture photometry on the sky noise images with the aperture correction to obtain the total flux f i . The steps 1)-3) are repeated for 100 times. We make a histogram of f i , and fit the histogram with a Gaussian profile. We regard the standard deviation divided by the median magnification factor of the subsample as the uncertainty of the total flux. To investigate the typical stellar mass M for a given M UV , we conduct the SED fitting method. We use the BEAGLE tool (Chevallard & Charlot 2016), which is based on a recent version of the stellar population models of Bruzual & Charlot (2003) and the photoionization models of Gutkin et al. (2016) that are computed with CLOUDY (Ferland et al. 2013). The intergalactic medium (IGM) absorptions follow the models of Inoue et al. (2014). We adopt the Calzetti et al. (1994) law to the models for dust attenuation.

SED Modeling

There are six free parameters in the fitting: i) redshift z of the galaxy, ii) galaxy age t age , iii) galaxy-wide ionization parameter U ion , iv) total mass of the formed stars M , v) stellar metallicity Z, and vi) V 606 -band dust attenuation optical depth τ V . We assume uniform prior probability distribution functions (PDFs) in the range of

z median phot -0.1 ≤ z ≤ z median phot + 0.1, (3)
1 Myr ≤ t age ≤ t univ (z = z median phot )t univ (z = 15), ( 4)

-3 ≤ log U ion ≤ -1, (5) 6 ≤ log(M /M ) ≤ 11, (6) -2 ≤ (log Z/Z ) ≤ 0.2, and (7)
0 ≤ log τ V ≤ 2 (8)
for the parameters i), ii), iii), iv), v), and vi), respectively. Here, z median phot is the median of the photometric redshifts z phot of the subsample whose values are taken from [START_REF] Kawamata | [END_REF]. We omit z phot values that are lower than 4 to derive z median phot in the same manner as the estimations of the median photometric redshifts [START_REF] Kawamata | [END_REF]. The age of the universe at z is represented by t univ (z). We assume a constant star formation history. The interstellar medium metallicity is assumed to equal the stellar metallicity. The dust-tometal ratio is fixed to 0.

L(Θ) = - 1 2 k f k -fk (Θ) σ k 2 . ( 10
)
Here f k , fk (Θ), and σ k are the observed flux, the flux predicted by the parameter set Θ, and the flux uncertainty, respectively. The subscript k runs over all of the bands but the i 814 , Y 105 , and J 125 band for z ∼ 6 -7, 8, and 9 subsamples, respectively. We do not include the i 814 , Y 105 , and J 125 band for z ∼ 6 -7, 8, and 9 subsamples, respectively. This is because we need to avoid the broadband photometry contaminated by unknown Lyα emission and IGM absorption effects. The value of σ k is defined as

σ k ≡ (σ obs k ) 2 + (σ 0 f k ) 2 , ( 11
)
where σ obs k is the observational uncertainty estimated in Section 3.1, and σ 0 is the relative systematical uncertainty, e.g., errors in background subtraction, flux calibration, and model predictions (Brammer et al. 2008;Dahlen et al. 2013;Acquaviva et al. 2015;Chevallard & Charlot 2016). We define σ 0 = 0.04 (0.05) for the HST and K s (IRAC) bands, applying the values recommended in the user manual of BEAGLE (Chevallard & Charlot 2016). If the S/N ratio of the stacked image in the band k is less than 2, we substitute f k values with 0. The posterior PDFs are efficiently sampled by the Nested Sampling algorithm implemented in the MULTINEST tool (Feroz & Hobson 2008;Feroz et al. 2009). We present the best-fit SEDs with the data photometries for z ∼ 6 -7, 8, and 9 subsamples in Figure 4, 5, and 6, respectively. An example of the posterior PDFs of the parameters is shown in Figure 7. We fit a linear function

Stellar Mass to UV Luminosity Relations

log(M /M ) = a 0 + a 1 (M UV + 19.5), (12)
where the intercept a 0 and the slope a 1 are set as free parameters, to our data in the magnitude range of -21.0 ≤ M UV ≤ -16.0. At z ∼ 8 and 9, we fix a 1 to the best-fit value of a 1 of z ∼ 6 -7, because fitting is unstable due to the lack of data points. The black solid lines in Figure 8 show the best-fit relations whose parameters are listed in Table 2. We find that the best-fit value of a 1 at z ∼ 6 -7 is -0.47 ± 0.09 that is similar to those obtained in the previous studies within the 1σ uncertainties. This value of a 1 is consistent with -0.4 that is the slope for the constant M /L UV case, indicating that the average M /L UV of dropouts at z ∼ 6 -7 is nearly constant at M ∼ 10 6 -10 9 M . The best-fit values of intercept a 0 are 8.67 +0.11 -0.11 , 8.59 +0.18 -0.19 , and 8.76 +0.13 -0.15 at z ∼ 6 -7, 8, and 9, respectively, showing no evolution beyond the errors. No evolution of a 0 indicates that stellar populations of the dropouts do not significantly change at z ∼ 6 -9 in the magnitude range.

In Figure 8, the blue dashed and dotted lines denote the M -M UV relations (corresponding to M /L UV ratios) that are assumed by Bouwens et al. (2017) for their M estimations from their M UV measurements for z ∼ 6 -8 dropouts. Bouwens et al. (2017)

respectively. Figure 8 indicates that our best-fit M -M UV relations (the black solid lines) are comparable to the model of 100-Myr star-formation duration time (the blue dotted lines). However, our best-fit M -M UV relations fall above the model of 10-Myr star-formation duration time (the blue dashed lines). By the comparison with our best-fit parameters at z ∼ 6 -7, we find that the parameter set for the blue dashed line (a 0 = 8.0 and a 1 = -0.4) is ruled out at the 5σ confidence level.

The parameter set of (a 0 , a 1 ) = (8.0, -0.4) provides the M -M UV relation 2 -7 times lower than the one of our best-fit parameters at z ∼ 6 -9 in M .

Galaxy Stellar Mass Functions

To derive the GSMFs, we apply the best-fit M -M UV relations (Section 3.3) to the LFs of Ishigaki et al. (2018). We show our GSMFs in Figure 9 with black dots. One can find that our GSMFs are broadly consistent with those of previous studies in the range of M ∼ 10 7 -10 9 M at z ∼ 6 -7. In the low-mass end, our GSMFs reach M ∼ 10 6 M , which is lower than those of the previous studies by 1 dex.

We parametrize our GSMFs with the Schechter (1976) function,

Ψ(M ) = ln(10)Ψ * M M * α+1 exp - M M * . ( 15
)
The characteristic stellar mass M * , the low-mass end slope α, and the normalization Ψ * are free parameters whose best-fit values are obtained by the least-squares fitting method. To improve the statistical accuracy, we fit Schechter functions simultaneously to our GSMFs and the GSMFs derived by Song et al. (2016) who do not use the HFF data but the CANDELS/GOODS (Grogin et al. 2011;Koekemoer et al. 2011) and the HUDF (Beckwith et al. 2006;Illingworth et al. 2013) data. We combine our GSMFs at z ∼ 6 -7 and the GSMFs of Song et al. (2016) at z ∼ 6, because the median value of z phot of our dropouts is z = 6.3 that is regarded as z ∼ 6. Table 3 summarizes the best-fit Schechter parameters at each redshift.

Galaxy Stellar Mass Density

We derive the galaxy stellar mass densities (GSMDs) at z ∼ 6 -9. We integrate our best-fit Schechter functions of the GSMFs over the stellar mass range of 10 8 ≤ M /M ≤ 10 13 that is adopted in previous studies. Table 4 lists our GSMDs at z ∼ 6 -9. The uncertainties of our GSMDs include the statistical errors and the cosmic variance uncertainties, the latter of which are added in quadrature with the values of 0.3, 0.4 and 0.5 dex at z ∼ 6 -7, 8 and 9, respectively (Robertson et al. 2014;Ishigaki et al. 2015;Bhatawdekar et al. 2018).

We show our GSMDs in Figure 10 together with the results obtained in the previous studies. Our GSMDs are consistent with those of the previous studies within 1σ uncertainties.

If the GSMDs and the star-formation rate densities (SFRDs) are accurately measured, the time integral of the SFRDs should be consistent with the GSMDs. In Figure 10, we show the evolution of the GSMDs that is derived as the time integration of the SFRDs obtained in Madau & Dickinson (2014) with the black solid curve. The SFRDs of Madau & Dickinson (2014) declines smoothly by (1 + z) -3.9 at z > 8. Similar trends of the SFRD evolution are found in the studies of Finkelstein et al. (2015), McLeod et al. (2016) and Bhatawdekar et al. (2018; see also Ellis et al. 2013 andMadau 2018). However, Oesch et al. (2018) claim that the SFRDs evolve more rapidly by (1 + z) -10.9 at z > 8 that is shown with the blue solid curve in Figure 10. There are two scenarios of the SFRD evolutionary trends, the smooth evolution of (1+z) -3.9 and the rapid evolution of (1 + z) -10.9 . We calculate χ 2 values with our GSMDs for the two functions of the smooth evolution and the rapid evolution. The χ 2 value is the smaller for the smooth evolution than for the rapid evolution by a factor of ≈ 2.8, implying that our results support the smooth evolution of the SFRDs. The smooth evolution of the SFRDs suggests that the star-formation rate per dark-matter mass increases at z > 8 (Harikane et al. 2018;Oesch et al. 2018).

Size-Mass Relations

We compare the distribution of the local GCs and our dropouts in the parameter space of R e and M . The GCs in the Milky Way today have R e values up to ≈ 40 physical pc (e.g. Pal 5 and Pal 14) and M up to ∼ 10 7 M (e.g. NGC 5139) that are shown in the catalog compiled by Baumgardt & Hilker (2018). We thus regard the dropouts with R e ≤ 40 physical pc and M ≤ 10 7 M as early GC candidates.

We plot the R e -M values of our dropouts at z ∼ 6-7, 8, and 9 in Figure 11, 12, and 13, respectively. The values of R e are taken from [START_REF] Kawamata | [END_REF], who obtain the best-fit value by the fit of ellipsoidal Sérsic profiles to the galaxy profiles that are observed with the lensing effects. We measure M UV of the dropouts in the same manner as those for the stacked images. The values of M UV are converted to M with the best-fit M -M UV relations obtained in Section 3.3.

We compare the R e -M distributions with those of the stellar systems in the local universe (Norris et al. 2014 and the references therein). We find that two of our dropouts, HFF5C-4260-1364 and HFF5C-4039-1566,4 meet the criteria of the early GC candidates. The parameters of the candidates are listed in Table 5.

As stated in the introduction, Bouwens et al. (2017) study the R e -M distributions of the z ∼ 6-8 dropouts. However, Bouwens et al. (2017) use only the first four HFF cluster (A2744, M0416, M0717, and M1149) data that were available at the time of their publication. No studies of R e -M distribution for z 6 are conducted in the rest of two clusters, A1063 and A370. Our early GC candidates of HFF5C-4260-1364 and HFF5C-4039-1566 are placed in the A1063 field. These early GC candidates are identified for the first time.

In Bouwens et al. (2017), there are 18 z ∼ 6 -8 dropouts meeting the criteria of the early GC candidates with their stellar masses obtained with the assumption of the 100-Myr star-formation duration time (Equation 0 14). One out of the 18 dropouts, M0416I-6115434445, is also selected in our study that refers to the dropout as HFF2C-1156-3446. This dropout is also reported as a compact star-forming region GC1 in Vanzella et al. (2017). The stellar masses and effective radii obtained in our study, Bouwens et al. (2017) The previous studies suggest that the GCs in the Milky Way form almost instantaneously at z 5, arguing that the stars inside the GCs have high number densities and tight isochrones on the HR diagram (e.g., Vandenberg et al. 1996;Forbes & Bridges 2010). The dropouts that we select are thus candidates of star clusters that should be a part or a dominant component of high-redshift low-mass galaxies, some of which may be related to GCs today. Our calculations for observational feasibilities suggest that these candidates can be spectroscopically confirmed with the James Webb Space Telescope and the Extremely Large Telescopes such as the Thirty Meter Telescope.

SUMMARY

We present stellar populations of low-mass dropout galaxies at z ∼ 6 -9 down to a stellar mass of M ∼ 10 6 M identified with the full data sets of the HFF. The stellar populations are studied with the HST and Spitzer deep optical to NIR images and the cluster-lensing magnification effects. The major results of our study are summarized below.

1. We investigate the stellar populations of the dropout galaxies with the optical/NIR photometric measurements and the BEAGLE SED modeling tool. We estimate M as a function of the UV magnitude M UV at M ∼ 10 6 -10 9 M . We find that our best-estimate M /L UV is consistent with the constant M /L UV within the 1σ uncertainties.

Our best-estimate M /L UV function is comparable to a model of star-formation duration time of 100 Myr, but 2 -7 times higher than the one of 10 Myr assumed in a previous study at the 5σ level.

2. We apply the best-estimate M -M UV relation (corresponding to M /L UV ) to the LFs of Ishigaki et al. (2018) to derive the GSMFs at z ∼ 6 -9. The GSMFs agree with those obtained by previous studies broadly at M 10 8 M , extending to the low-mass limit of ∼ 10 6 M . We reach this lowmass limit, exploiting the deep HFF full data sets and the gravitational-lensing magnification effects.

3. Integrating the GSMFs, we find that the GSMDs ρ smoothly increases from log(ρ /[M Mpc -3 ]) = 5.91 +0.75 -0.65 at z ∼ 9 to 6.21 +0.39 -0.37 at z ∼ 6 -7. This trend is consistent with the one estimated from time integration of the SFRD evolution compiled by Madau & Dickinson (2014).

4. The stellar masses of the dropouts are presented as a function of effective radius R e with our best-fit M -M UV relations. We find that the two dropouts have low stellar masses (M ≤ 10 7 M ) and compact morphologies (R e ≤ 40 physical pc), which are comparable with those of GCs in the Milky Way today. These dropouts are candidates of star clusters that should be a part or a dominant component of high-redshift low-mass galaxies. These dropouts may be related to the present GCs, given the fact that the old stellar populations are found in GCs of the Milky Way today. show the distributions of the local elliptical/S0 galaxies (E/S0s; black circles), dwarf elliptical/S0 galaxies (dEs/dS0s; gray triangles), dwarf spheroids (dSphs; light-green squares), nuclear star clusters (NSCs; orange crosses), young massive clusters (YMCs; green down-pointing triangles), and globular clusters/ultra-compact dwarfs (GCs/UCDs; blue hexagons) obtained by Norris et al. (2014) and the references therein. The cyan shade represents the region where the local GCs are located. We define the dropouts in this region as the early GC candidates. The right vertical axis shows Re in the angular scale at z = 6.

M-Figure 1 .-Figure 2 .-Figure 3 .

 123 Figure 1. Stacked images of the z ∼ 6 -7 subsamples for each band. The image size is 4 × 4 .

  z ∼ 6 -7, MUV = -21.3 MUV ≤ -21.0 4 -20.6 -21.0 ≤ MUV ≤ -20.5 7 -20.2 -20.5 ≤ MUV ≤ -20.0 19 -19.7 -20.0 ≤ MUV ≤ -19.5 43 -19.2 -19.5 ≤ MUV ≤ -19.0 58 -18.8 -19.0 ≤ MUV ≤ -18.5 64 -18.3 -18.5 ≤ MUV ≤ -18.0 30 -17.7 -18.0 ≤ MUV ≤ -17.5 18 -17.3 -17.5 ≤ MUV ≤ -17.0 15 -16.6 -17.0 ≤ MUV ≤ -16.0 10 -15.1 -16.0 ≤ MUV 5 z ∼ 8, MUV = -20.4 MUV ≤ -20.0 13 -19.7 -20.0 ≤ MUV ≤ -19.5 10 -19.3 -19.5 ≤ MUV ≤ -19.25 10 -19.1 -19.25 ≤ MUV ≤ -19.0 11 -18.7 -19.0 ≤ MUV 10 z ∼ 9, MUV = -20.3 MUV ≤ -19.9 8 -19.6 -19.9 ≤ MUV ≤ -19.5 9 -19.3 -19.5 ≤ MUV ≤ -19.1 9 -18.7 -19.1 ≤ MUV 10 Note-Columns: (1) subsample name that indicates the redshift and the median value of MUV; (2) threshold of the subsample; (3) number of the dropouts in the subsample.

Figure 8

 8 Figure 8 shows the relations between the stellar mass M and the UV magnitude M UV . Here, M and M UV represent the median value of the marginal posterior PDF and the median value of the UV magnitudes of the dropouts in the subsample, respectively. The vertical error bars represent 68 % confidence intervals, while the horizontal bars show the minimum and maximum values of M UV of the dropouts in the subsamples. Figure 8 indicates that our results are broadly consistent with the previous results.We fit a linear function

Figure 4 .

 4 Figure 4. Intrinsic SEDs the z ∼ 6 -7 subsamples. The subsample names are denoted at the upper-left corner in each panel. The red filled squares show the photometric data obtained in Section 3.1. The red open squares at i814band are the same as the red filled squares, but for the photometry that are not included in the likelihood calculation (see text). The horizontal and vertical error bars represent the wavelength range of the filters and the 1σ uncertainties, respectively. The down arrows indicate 2σ upper limits for bands with S/N < 2. The black lines represent the bestfit SEDs, while the blue crosses show the bandpass-averaged magnitudes predicted from the best-fit SEDs. In the bottom panels, the color shades denote the normalized filter throughputs of B435, V606, i814 (cyan), Y105, J125, JH140, H160 (green), Ks (orange), ch1, and ch2 (red).

Figure 5 .Figure 6 .

 56 Figure 5. Same as Figure 4, but for the z ∼ 8 subsamples.

Figure 7 .

 7 Figure 7. Posterior PDFs of the parameters for the z ∼ 6 -7, MUV = -21.3 subsample. The diagonal panels show the marginal distributions of z, log(M /M ), log(tage), log(Z/Z ), and AV. Note that AV is derived from τV. The red vertical lines and the shades represent the median values and the 1σ intervals, respectively. The off-diagonal panels represent the joint distributions for the sets of two parameters. The inner, middle, and outer contours correspond to 1σ, 2σ, and 3σ intervals, respectively.

Figure 8 .

 8 Figure 8. Distribution of M versus MUV of the z ∼ 6 -7 (left), 8 (middle), and 9 (right) subsamples. The red filled circles represent our results. The black lines and the gray shades denote the best-fit relations and the 1σ uncertainties, respectively. The blue dashed (dotted) lines show the M -M relation of the representative model of the 10 (100) Myr star-formation duration time that is assumed in Bouwens et al. (2017). The other data points are taken from the previous studies, Stark et al. (2013, orange triangles), Duncan et al. (2014, cyan diamonds), Song et al. (2016, purple pentagons), and Bhatawdekar et al. (2018, green boxes). In the left panel, the filled and open symbols show the data at z ∼ 6 and 7, respectively. All of the IMFs are converted to a Chabrier (2003) IMF.

Figure 9 .

 9 Figure 9. GSMFs at z ∼ 6 -7 (left), 8 (middle), and 9 (right). The red filled circles show our GSMFs. The red open circles are the same as the red filled circles, but for the GSMFs obtained with the extrapolated best-fit M -MUV relations. We do not include the data of the red open circles for our Schechter function fitting. The down arrows denote the upper limits of the GSMFs. The other data points are taken from the previous studies, Bhatawdekar et al. (2018, green boxes), Song et al. (2016, orange pentagons), Grazian et al. (2015, black crosses), Duncan et al. (2014, cyan diamonds), and González et al. (2011, magenta triangles). In the left panel, the filled and open symbols show the GSMFs at z ∼ 6 and 7, respectively. The black curves represent the best-fit Schechter functions obtained with the measurements of our study and Song et al. (2016). The stellar masses taken from the previous studies are converted to those estimated with the Chabrier (2003) IMF.

  : K19 = this work, D14 = Duncan et al. (2014), G15 = Grazian et al. (2015), S16 = Song et al. (2016), and B18 = Bhatawdekar et al. (2018). Here we show the 'point source' results fromBhatawdekar et al. (2018).

Figure 10 .

 10 Figure 10. GSMDs as a function of redshift. The red circles show our GSMDs. Here the error bars include the statistical errors and cosmic variance uncertainties. Our GSMD at z ∼ 6 -7 is plotted at z = 6.5. The other data points are taken from the previous studies, Bhatawdekar et al. (2018, green diamonds), Song et al. (2016, purple boxes), Grazian et al. (2015, cyan pentagons), Oesch et al. (2014, a blue diamond), Duncan et al. (2014, orange stars), Stark et al. (2013, cyan hexagons), González et al. (2011, gray hexagons), Mortlock et al. (2011, magenta diamonds), and Elsner et al. (2008, navy triangles). The data points are slightly shifted along the redshift axis if overlapped with other points. The black and blue solid curves indicate the GSMD evolution that is derived as the time integration of the SFRD evolution presented in Madau & Dickinson (2014) and Oesch et al. (2018), respectively. The stellar masses taken from the previous studies are converted to those estimated with the Chabrier (2003) IMF.

Figure 11 .

 11 Figure11. Effective radius on the source plane, as a function of M . Our i-dropouts are presented with the large circles whose colors indicate the magnification factor values that are defined with the color bar at the bottom right. The other symbols show the distributions of the local elliptical/S0 galaxies (E/S0s; black circles), dwarf elliptical/S0 galaxies (dEs/dS0s; gray triangles), dwarf spheroids (dSphs; light-green squares), nuclear star clusters (NSCs; orange crosses), young massive clusters (YMCs; green down-pointing triangles), and globular clusters/ultra-compact dwarfs (GCs/UCDs; blue hexagons) obtained byNorris et al. (2014) and the references therein. The cyan shade represents the region where the local GCs are located. We define the dropouts in this region as the early GC candidates. The right vertical axis shows Re in the angular scale at z = 6.

Figure 12 .Figure 13 .

 1213 Figure 12. Same as Figure 11, but for Y -dropouts. The right vertical axis shows Re in the angular scale at z = 8.

  3 (e.g., De Vis et al. 2017).

The posterior PDF P (Θ|D, H) of a given parameter set Θ = {z, t age , U ion , M , Z, τ V } of a model H is calculated based on the Bayes' theorem P (Θ|D, H) = P (Θ|H)P (D|Θ, H) dΘP (Θ|H)P (D|Θ, H) (9) (e.g., Jeffreys 1961), where D is the data set, i.e., the fluxes in the B 435 , V 606 , i 814 , Y 105 , J 125 , JH 140 , H 160 , K s , ch1, and ch2 bands. The likelihood function of Θ, L(Θ) ≡ P (D|Θ, H), is defined via ln

Table 2 .

 2 Best-fit parameters of the M -MUV relations.

	z	a0	a1
	6 -7 8.67 +0.11 -0.11	-0.47 +0.09 -0.09
	8	8.59 +0.18 -0.19	-0.47 (fixed)
	9	8.76 +0.13 -0.15	-0.47 (fixed)

Table 3 .

 3 Best-fit parameters of the Schechter function of the GSMFs.

	Ref.	M	α	Ψ *
		[log M ]		[10 -5 dex -1 Mpc -3 ]
	z ∼ 6
	K19	9.58 +0.23 -0.15	-1.85 +0.07 -0.07	18.2 +12.7 -9.3
	D14	10.87 +1.13 -1.06	-2.00 +0.57 -0.40	1.4 +41.1 -1.4
	G15 10.49 ± 0.32 -1.55 ± 0.19	6.19 +13.50 -4.57
	S16	10.72 +0.29 -0.30	-1.91 +0.09 -0.09	1.35 +1.66 -0.75
	B18	10.29 +0.48 -0.56	-1.93 +0.05 -0.07	5.63 +7.12 -3.23
	z ∼ 7
	D14 10.51 (fixed)	-1.89 +1.39 -0.61	3.6 +30.1 -3.5
	G15 10.69 ± 1.58 -1.88 ± 0.36	0.57 +59.68 -0.57
	S16	10.78 +0.29 -0.28	-1.95 +0.18 -0.18	0.53 +1.10 -0.38
	B18	10.25 +0.45 -0.49	-1.95 +0.06 -0.07	2.82 +5.16 -1.88
	z ∼ 8
	K19	8.93 +0.23 -0.15	-1.52 +0.27 -0.26	45.7 +33.7 -25.7
	S16	10.72 +0.29 -0.29	-2.25 +0.72 -0.35	0.035 +0.246 -0.030
	B18	10.48 +0.90 -0.75	-2.25 +0.23 -0.29	0.089 +0.36 -0.074
	z ∼ 9
	K19	9.04 +0.47 -0.18	-1.55 +0.29 -0.30	60.3 +37.5 -46.8
	B18	10.45 +0.80 -0.85	-2.33 +0.30 -0.39	0.06 +0.66

Table 4 .

 4 Galaxy stellar mass densities.

	z	log ρ [M Mpc -3 ]
	6 -7	6.21 +0.39 -0.37
	8	5.64 +0.52 -0.51
	9	5.91 +0.75 -0.65

Table 5 .

 5 Catalog of the early GC candidates.

	ID	R.A. (J2000) decl. (J2000) MUV [mag] M [10 6 M ] Re [physical pc]	µ
	HFF5C-4260-1364	342.177526	-44.526786 -15.7 ± 0.2	7.7 ± 7.1	40 +30 -20	19.12
	HFF5C-4039-1566	342.168321	-44.532412 -14.2 ± 0.1	1.6 ± 2.0	40 +30 -20	77.69

http://cosmos.phy.tufts.edu/ ∼ danilo/HFF/Download.html

https://archive.stsci.edu/prepds/frontier/lensmodels/

These IDs are identical to the ones by[START_REF] Kawamata | [END_REF] andIshigaki et al. (2018).

We are grateful to Michiko Fujii, Seiji Fujimoto, Ryo Higuchi, Ryohei Itoh, Nobunari Kashikawa, Shiro Mukae, and Haibin Zhang for useful comments and discussions. We particularly thank Masafumi Ishigaki and Ryota Kawamata for their sample catalogs and mass models. GB acknowledges financial support through PAPIIT projects IG100319 from DGAPA-UNAM. This work is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (15H02064, 17H01110, 17H01114, 15H05892, and 18K03693) Grant-in-Aid for Scientific Research (A) through Japan Society for the Promotion of Science. This work is based on data and catalog products from HFF-DeepSpace, funded by the National Science Foundation and Space Telescope Science Institute (operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555).


B Desplanques
email: desplanq@lpsc.in2p3.fr

Y B Dong
email: dongyb@mail.ihep.ac.cn

Form factors in relativistic quantum mechanic: constraints from space-time translations

The comparison of form factors calculated from a single-particle current in different relativistic quantum mechanic approaches evidences tremendous discrepancies. The role of constraints coming from space-time translations is considered here with this respect. It is known that invariance under these translations implies the energy-momentum conservation relation that is usually assumed to hold globally. Transformations of the current under these translations, which lead to this result, also imply constraints that have been ignored so far in relativistic quantum mechanic approaches. An implementation of these constraints is discussed in the case of a model with two scalar constituents. It amounts to incorporate selected two-body currents to all orders in the interaction. Discrepancies for form factors in different approaches can thus be removed, contributing to restore the equivalence of different approaches. Results for the standard front-form approach (q + = 0) are found to fulfill the constraints and are therefore unchanged. The relation with results from a dispersion-relation approach is also made.

Introduction

Examination of form factors calculated in different forms of relativistic quantum mechanic (RQM) [1,2], with the same solution of a mass operator, shows tremendous differences depending on the approach that is used [3,4]. On the other hand, fitting the wave function to experiment shows that it depends strongly on the approach [5,6]. At first sight, one cannot therefore decide whether some discrepancy is due to the chosen implementation of relativity or to the underlying dynamics. It is however believed that results should not depend on the choice of the hypersurface underlying some form of relativity [7] and that their equivalence requires the introduction of many-body currents.

Among properties expected from Poincaré covariance, invariance of form factors under rotations or boosts can easily be checked by applying these transformations to the system under consideration. In absence of a similar test, the role of space-time translations, which are also part of the Poincaré group, is much less known beyond the energy-momentum conservation that is, of course, assumed to hold. In this contribution, we examine this role and show that accounting for the constraints they imply can remove large discrepancies between results obtained in different approaches with the simplest one-body current.

In the following, we successively consider the space-time translations and constraints they imply (sect. 2), the implementation of these constraints (sect. 3), the illustration of their role for a pion-like system (sect. 4) and, finally, the relation to a dispersion-relation approach that could be considered as the convergence point of all approaches based on the same solution of a mass operator (sect. 5). Due to lack of space, details had to be skipped. They could be found in ref. [8].

Constraints from Poincaré space-time translation invariance

Under space-time translations, a vector or a scalar current transforms as follows:

e iP •a J ν (x) (or S(x)) e -iP •a = J ν (x + a) (or S(x + a)), (1)
where P µ represents the total momentum operator. Considering the matrix element of this equality between eigenstates of P µ for a = -x, one gets:

< i |J ν (x) (or S(x))| f >= e i (P i -P f )•x < i |J ν (0) (or S(0))| f > . (2)
Together with a field carrying momentum q µ , one gets under the assumption of space-time translation invariance the current momentum-energy conservation, (P

f -P i ) µ = q µ .
For the purpose of calculating form factors, J ν (0) (or S(0)) is generally approximated by a one-body current. Equation (1) however implies further relations involving the commutator of P µ with the currents and the derivative of the current with respect to x [9]. Particularly interesting relations are the following double commutators:

P µ , P µ , J ν (x) = -∂ µ ∂ µ J ν (x), P µ , P µ , S(x) = -∂ µ ∂ µ S(x) . (3)
Considering the matrix element of these relations between eigenstates of P µ , and for the case of a one-body current, one should verify the relation:

< |q 2 J ν (0) (or S(0))| >=< |(p i -p f ) 2 J ν (0) (or S(0))| > . (4)
It is easily seen that this equation cannot be generally fulfilled, as most often q 2 = (p i -p f ) 2 in RQM approaches (see fig. 1 for a graphical representation). This implies that the assumption of a single-particle current is inconsistent with properties from space-time translations and that the current, J ν (0) (or S(0)), besides a one-body component, should also contain many-body components which, until now, have been ignored.

E , P

Implementation of the constraints

It is expected that many-body currents at all orders in the interaction are required to fulfill constraints from space-time translations. As considering these currents explicitly is excluded, we consider them indirectly, by modifying wave functions and the current operator. The modification is suggested by examining expressions of form factors in various approaches, which show that the factor multiplying Q could be given by a factor varying from 1 in some cases to 2e k M in other ones [10]. The departure of this quantity to 1 represents an interaction term and is a signature of the hypersurface underlying the approach. To account for the constraints, we thus propose to multiply Q by a factor α and determine this one by requiring that the squared momentum transferred to the system, q 2 , be equal to the one for the constituents, denoted "(p i -p f )" 2 . The equation to be solved is typically given by:

q 2 = "[(P i -P f ) 2 + 2 (∆ i -∆ f ) (P i -P f ) • ξ + (∆ i -∆ f ) 2 ξ 2 ]" = α 2 q 2 -2 α "(∆ i -∆ f )" q • ξ + "(∆ i -∆ f ) 2 " ξ 2 , (5)
where ∆, which represents an interaction effect, also depends on α. Explicit expressions of α can be found in ref. [8] for different forms. Expressions for form factors, taking into account the effect of constraints motivated by space-time translation properties, can then be obtained [8]. By expanding these expressions in terms of ∆, the many-body character of corrections at all orders of the interaction could be checked.

A few points deserve to be noticed. Firstly, for the standard front-form where ξ 2 = 0, q • ξ (or q + ) = 0, the factor α is equal to 1 and, therefore, results for the form factors are unchanged. Secondly, a solution has been found but we do not exclude that the implementation can be done differently, what would be desirable to fulfill the infinite set of constraints involving the multiple commutators of P µ with currents. Thirdly, the choice of the charge current is constrained for some part. We used for our purpose a current inspired from results for the simplest Feynman triangle diagram. This allows one to get Lorentz invariant results for form factors. To illustrate effects of the restoration of properties related to space-time translations, we consider a system of scalar particles interacting through the exchange of massless particles (Wick-Cutkosky model [11,12]). For this system, that may be used as a test case [13], form factors can be calculated exactly. There are two of them for the ground state, a charge and a scalar one, which provide us with some "experiment". The RQM calculations are based on the solution of a mass operator where the interaction is chosen so that to reproduce approximately the degeneracy pattern of the exact spectrum while keeping the high-momentum power law unchanged. The strength is fitted to reproduce the energy of the ground state used in the calculations. The mass of the system and of the constituents are those currently used for the pion, M = 0.14 GeV and m = 0.3 GeV.

Results are presented in figs. 2, 3 for the charge and scalar form factors. In each case, there are two panels describing the low-and high-Q 2 behavior, respectively sensitive to the square radius and to the asymptotic behavior. The different curves correspond to uncorrected Breit-frame form factors for the front form (q + = 0; F.F. (perp.)), a front form with a parallel configuration (F.F. (parallel)), the usual instant form (I.F.) [14], some instant form with the symmetry properties of the point form ("P.F.") [15], a point form inspired from the Dirac one (D.P.F.) [16], and corrected results which, all, coincide with the standard front-form results (q + = 0) that are unchanged. Anticipating on next section, we also show results of a dispersion-relation approach (D.R.), which coincide with the standard front-form ones. Examination of results for the charge form factor (fig. 2) shows that tremendous discrepancies are removed by accounting for constraints from space-time translation properties at both low and high Q 2 . The paradox of a charge radius tending to infinity while the mass of the system goes to zero (or the interaction is increased) disappears. There remain a slight discrepancy with "experiment". It can be ascribed for a part to the description of the mass operator. It is not clear whether genuine two-body currents are needed. Examination of results for the scalar form factor (fig. 3) shows similar features as far as corrections related to space-time translation properties are concerned. In contrast to the charge form factor, genuine two-body currents are needed here to explain the form factor at low Q 2 . It is reminded that, contrary to the charge form factor, the scalar one at Q 2 = 0 is not protected by some conservation law.

Relation to a dispersion-relation approach

A dispersion-relation has been proposed to calculate form factors [17,18,19]. The expressions for the charge and scalar form factors read:

F 1 (Q 2 ) = 1 N ds d s i -s f Q (2s + Q 2 ) θ s i s f D -m 2 D 3/2 φ(s f ) φ(s i ) , F 0 (Q 2 ) = 1 N ds d s i -s f Q θ s i s f D -m 2 2 D 1/2 φ(s f ) φ(s i ) ,
where

N = ds s -4 m 2 s φ 2 (s) , s = s i +s f 2 , D = 4s+Q 2 + (s i -s f ) 2 Q 2 . (6)
This approach was presented as an instant-form one [18]. In absence of reference to a particular direction, it could be thought as a point-form one. Moreover, one of the authors (B.D.) found that the above expressions could be obtained from the standard front-form ones (q + = 0) by an appropriate change of variables. A similar result was obtained independently by Melikhov [19] for the pion case. Examining the approach, it is noticed that it is based on the free scattering amplitude in an external field. The interaction effects that are here or there in RQM approaches are therefore absent. The variables s i = (p i + p) 2 , s f = (p f + p) 2 and (p i -p f ) 2 refer to on-mass shell constituents.

The functions φ(s) = φ(k 2 = s 4 -m 2 ) can therefore be identified to the solution of a mass operator used in RQM approaches. The relation q 2 = (p i -p f ) 2 , which ensures that the square momentum transferred to the system be equal to the one transferred to the constituents, fulfills the constraints expected from space-time translation properties. These three features suggest that the RQM results for form factors should converge to the dispersion-relation ones, once the above constraints are accounted for.

As the relations of corrected RQM form factors to the dispersion-relation ones are not well known, we give here expressions evidencing the relation for the charge form factor. The explicit form of the change of variables can be found in ref. [8]. The first case refers to the standard front form approach (q + = 0), which is of particular interest for this light-cone meeting. Using the k ⊥ and x variables currently employed in this domain, the relation reads:

F 1 (Q 2 ) = 2 πN d 2 p ⊥ dx 2x(1-x) φ( k f 2 ) φ( k i 2 ) = 2 πN ds d( s i -s f Q ) 1 2 √ D φ(s f ) φ(s i ) × dx (2 s+Q 2 )/D -(x-d) s i s f D -m 2 f -(x-d) 2 = 1 N ds d( s i -s f Q ) (2 s + Q 2 ) θ(• • •) D 3/2 φ(s f ) φ(s i ) . (7)
The second case refers to an arbitrary front orientation λ µ with finite λ 2 . The relation reads:

F 1 (Q 2 ) = 16π 2 N(2π) 3 d p e p "(2p+p i +p f )"•λ 2"(p i +p f )"•λ " φ( k f 2 ) φ( k i 2 )" = 2 πN ds d( s i -s f Q ) θ(• • •) 4 D 3/2 φ(s f ) φ(s i )× d(p• λ) (2s+Q 2 ) -(p• λ-d)g s i s f D -m 2 f -(p• λ-d) 2 . (8)
In this case, it is noticed that the dependence on λ µ and P µ = (P i + P f ) µ /2, which characterizes quantities d, f and g [8], disappears in the integration over p • λ, allowing one to get Lorentz-invariant results. In both cases, the integral is reduced from a threeto a two-dimensional one.

Conclusion

In this contribution, we mainly considered properties related to Poincaré space-time translations in RQM approaches for the calculation of form factors. The current practice is to use these properties to factorize the x dependence of the current, allowing one to get the usual energy-momentum conservation relation, and to assume that the current at x = 0 is a one-body current. At first sight, there is no relation of the 4-momentum transferred to the system and to the constituents as in a field-theory approach. Considering further relations implied by the transformation of the current under space-time translations, we showed that they could not be satisfied with a current reduced to a one-body component, pointing to the necessary presence of many-body components. These ones allow us to account in RQM approaches for the equality of the momentum transferred to the whole system and to the constituents, which characterizes a field-theory approach. We described a method to implement the above many-body contributions in the case of a scalar-particle model. When this is done, it is found that discrepancies between different RQM approaches for calculating form factors can be removed, showing that the role of space-time translations extends beyond the standard energy-momentum conservation.

It is also found that these results could coincide with those of a dispersion-relation approach. Altogether, all aspects of the Poincaré group (rotations, boosts and space-time translations) are essential in getting reliable estimates of form factors as far as the implementation of relativity is concerned. Possible discrepancies with experiment can thus be more likely ascribed to the underlying dynamics.

Present results can be extended without much difficulty to elastic as well as inelastic form factors of 0-spin systems consisting of 1/2-spin constituents with unequal masses. Extension to non-zero spin systems could require more elaboration.

One of the authors (B.D.) is very grateful to V.A. Karmanov and W. Polyzou for an initiation to implementation of relativity in field theory and relativistic quantum mechanic respectively. It is likely that the complementary knowledge of these two approaches was essential in obtaining results presented here. This work is partly supported by the National Sciences Foundations of China under grant No. 10775148.

Figure 1 :

 1 Figure 1: Representation of the interaction with an external probe

Figure 2 :

 2 Figure 2: Charge form factors. See text for the curves. The "experiment" is represented by diamonds.

Figure 3 :

 3 Figure 3: Scalar form factors. See text for the curves. The "experiment" is represented by diamonds.


Maxime Martinez

Olivier Giraud

Denis Ullmo

Juliette Billy

David Guéry-Odelin

Bertrand Georgeot

Gabriel Lemarié

Chaos-assisted long-range tunneling for quantum simulation

come

Introduction.-In recent years there has been considerable interest in the quantum simulation of more and more complex problems of solid state physics [1][2][3]. In this context, lattice-based quantum simulation has become a key technique to mimic the periodicity of a crystal structure. In such systems, dynamics is governed by two different types of processes: hopping between sites mediated by tunneling effect and interaction between particles. While there exists several ways to implement long-range interactions [4][5][6][7], long-range hoppings have been up to now very challenging to simulate. These long-range hoppings however, have aroused great theoretical interest in condensed matter, as they are associated with important problems such as glassy physics [8], many-body localization [9] or quantum multifractality [10]. In this study we show that such long-range hoppings can be engineered in driven lattices in a moderate regime of modulation.

Temporal driving techniques are widely used in quantum simulation [11], as systems with fast driving can exhibit new topological effects [12][13][14][15][16] and systems with strong driving can mimic disorder [START_REF] Casati | Stochastic Behavior in Classical and Quantum Hamiltonian Systems[END_REF][18][19][20][21][22]. In the intermediate regime of modulation we focus on, driven lattices have a classical dynamics which is neither fully chaotic (corresponding to the strong driving case) nor regular (corresponding to the fast driving case), but, as most real-life dynamical systems, show coexistence of chaotic and regular zones in phase space. Our main result is based on the richness of the quantum tunneling effect in such systems: it is known to be chaos-assisted [23][24][25][26][27][28][29][30][31][32][33][34][START_REF] Keshavamurthy | Dynamical Tunneling -Theory and experiment[END_REF]. This phenomenon is well understood between two regular islands, where it translates into large resonances of the tunneling rate between the two islands when varying a parameter of the system. It has been observed in different experimental contexts, with electromagnetic waves [27,[START_REF] Dembowski | [END_REF][37][38][39][40][41], and with cold atoms [42][43][44][45].

In this paper, we address the generalization of this chaos-assisted tunneling mechanism to the case of a mixed lattice composed of regular islands embedded in a chaotic sea, obtained in a moderate regime of tem- poral driving. We demonstrate that remarkably such a dynamical quantum system can be mapped onto an effective tight-binding Hamiltonian with long-range hoppings ∝ 1/n, with n the distance between sites. Beyond the intrinsic interest of a new observable quantum chaos effect, our results open new engineering possibilities for latticebased quantum simulations as they are highly generic, accessible for state-of-the-art experiments and species independent (in the context of cold atoms).

V (x, t) (a) N λ (c) t n ∝ 1/n |0 |1 |n |N -2 |N -1 (b)
Model.-We consider an experimental situation similar to [45], that is to say a dilute condensate of cold atoms subjected to an optical lattice whose intensity is modulated periodically in time. Using dimensionless variables [START_REF]x = 2πX/d where P and X are the momentum and position along the optical lattice, M is the atomic mass, d is the lattice spacing, and the time t is measured in unit of the modulation angular frequency ω[END_REF], the dynamics is governed by the Hamiltonian

H(x, t) = p 2 2 -γ(1 + ε cos t) cos x. ( 1
)
Here γ is the dimensionless depth of the optical lattice and ε the amplitude of modulation. Importantly, the effective Planck constant heff = -i[x, p] = 2ω L /ω can be tuned experimentally (ω being the modulation angular frequency and ω L is related to the lattice characteristic energy E L = hω L ). This Hamiltonian has a dimensionless time period T = 2π and spatial period λ = 2π. Besides the choice of this model, our conclusions are valid for almost any modulation waveform (including phase modulation and kicked potentials).

Semiclassical picture.-The classical dynamics of this time-periodic system is best understood through a stroboscopic phase space, using the values of (x, p) at each period of the modulation t = jT , j integer. For ε = 0, the system is integrable. When ε increases, chaos develops close to the separatrix, forming a chaotic sea which surrounds regular islands of oscillating orbits centered around the stable points (x = 2nπ, p = 0, n an integer) of the potential wells, see Fig. 1. At ε = 0, with no chaotic sea, tunneling essentially occurs between adjacent wells, and the system can be described for deep optical lattices by an effective tight-binding Hamiltonian with nearest-neighbor hopping. Our main objective is to describe in a similar way the modulated system, a dynamical, spatially periodic lattice of N regular islands indexed by n ∈ [[0, N -1]], surrounded by a chaotic sea.

Adopting a stroboscopic point of view, the quantum dynamics is described by the evolution operator U F over one period of modulation. Each eigenstate |φ l of U F can be associated with a quasi-energy ε l , so that U F |φ l = exp(-iε l T /h eff ) |φ l . Equivalently the Hamiltonian H strob ≡ i(h eff /T ) log U F gives the exact same stroboscopic dynamics as U F and is associated with the same eigenstates |φ l with energies ε l .

In the semiclassical regime where heff < A, with A the area of a regular island, the quantum dynamics is strongly influenced by the structures of the classical phase space. Quantum eigenstates can be separated in two types [24,[START_REF] Berry | [END_REF]: regular (localized on top of regular orbits) or chaotic (spread over the chaotic sea), see Fig. 2 for a Husimi [48][49][50] phase-space representation.

The tunnel coupling between regular states is well understood in the case of N = 2 regular islands surrounded by a chaotic sea and is called chaos-assisted tunneling [23,24]. In the absence of a chaotic sea, tunneling involves only a doublet of symmetric and anti-symmetric states localized on the two symmetric islands. In the presence of a chaotic sea, the key property of chaosassisted tunneling is a 3-level mechanism with one of the regular states interacting resonantly with a chaotic state. This coupling translates in an energy shift of the involved regular state and thus of a strong variation of the original doublet energy splitting (which is nothing but the tunneling frequency). As a signature of the chaotic dynamics, it was shown [23,25] that these chaos-assisted resonances, observed for the first time in a quantum system only recently [45], occur quite erratically when varying a parameter of the system. Chaos-assisted tunneling involves a purely quantum transport (tunneling to the chaotic sea) and a classically allowed transport (diffusion in the chaotic sea). Thus in mixed lattices, long-range hopping can be expected because the chaotic sea connects all the regular islands across the lattice (see Fig. 1).

Effective Hamiltonian.-The existence of regular islands in the center of each cell (see Fig. 1) motivates the introduction of a set of regular states {|n reg } forming a lattice with one site per cell, whose exact construction is not crucial for our discussion (see [27] for a detailed discussion). For simplicity, we work in the regime heff < ∼ A such that there is only one regular state by island.

In clear contrast with regular lattices, where only neighboring sites are directly coupled by standard tunneling effect, there exists an indirect transfer of probability between distant sites of the modulated lattice due to the mutual coupling with additionnal chaotic states, strongly delocalized along the system. In the original scenario [23], the chaos-assisted tunneling mechanism between regular islands is associated with fast, but weak, probability oscillations between regular islands and the chaotic sea. This picture motivates to capture the physics of tunneling in our system through an effective Hamiltonian H eff , acting only in the regular subspace but generating the same dynamics as H strob in this subspace [51][52][53]. This effective Hamiltonian can be defined from (E -H eff ) -1 ≡ P reg (E -H strob ) -1 P reg , with P reg the projector onto the regular subspace spanned by the |n reg .

β0 Quasi-momentum β Quasi-energy ε(β) α |W | ∆β = 4|W | α 0 1 β±|βreg 2 -π/λ π/λ β ε(β) ≈ |βreg 1 √ 2 (|β ch + |βreg ) ≈ |β ch
In the effective picture, coupling with chaotic states translates in a shift of the energy of each regular Bloch state |β reg = 1 √ N n exp(iβλn) |n reg (with β an integer multiple of 2π/λN ). The resulting dressed regular band ε eff reg (β) then gives access to the effective tunneling coupling t eff n ≡ (m + n) reg |H eff |m reg through the Fourier transform in quasi-momentum

t eff n = 1 N β ε eff reg (β) exp(iβλn). ( 2
)
The simplest way to determine the effective spectrum is to start from the full exact spectrum (obtained numer-ically for instance) and to choose the N most relevant energies. The natural choice is to select energies associated with eigenstates with the largest projection on the regular subspace. In mixed lattices, this procedure gives systematic discontinuities in the effective band, coming from accidental degeneracies between a regular |β reg and a chaotic state |β ch . Close to such avoided crossings, the branch giving the effective regular energy suddenly changes, resulting in a sharp discontinuity of ε eff reg (β) (see Fig. 2). These sharp discontinuities cause, from property of the Fourier transform in Eq. ( 2), a long-range decay of the effective coupling term t eff n ∼ 1/n (see Fig. 4). The two main features of these resonances come from the mixed nature of the system. First, they are sharp because the local slope α = dε ch /dβ of the crossing state is large since ergodic chaotic states are sensitive to boundary conditions. Second, their heights 2|W | = 2| β ch |H|β reg | is larger than the regular band width 2| n reg |H|n + 1 reg | (nearest-neighbor hopping amplitudes for ε = 0, with direct tunneling between islands without chaotic sea). Numerical simulations.-To test the accuracy of this effective tight-binding picture, we compare the exact stroboscopic dynamics with the one given by the effective Hamiltonian, considering a wave packet initially localized on a single regular island of the modulated lattice. (see the Supplemental Material [54] for computational details). As concerns the exact dynamics, the initial condition was chosen to be a localised (Wannier) state of the undriven lattice (ε = 0), in the cell n 0 = (N -1)/2, N being odd. We also used the localized states |n reg to estimate the projection of the wavefunction on the chaotic layer through P ch ≡ 1 -P reg . The effective dynamics was studied by propagation of a state initially localized in the cell n 0 of a discrete lattice ruled by the effective Hamiltonian. The effective Hamiltonian was extracted from diagonalization of the Floquet operator on a single cell, for different values of β. In both types of simulations, we defined two observables: a local one pn which enables to probe the probability at each site, defined as pn ≡ |n n| in the effective system and pn ≡ (n+1)λ nλ |x x| dx in the exact system, (this choice was motivated to get n n = 1 in both systems) and a global one ∆n 2 = n (nn 0 ) 2 pn to estimate the spreading of the wave function.

We have simulated different system sizes up to N = 1079 with periodic boundary conditions and found a very good agreement between the two approaches (see Fig. 3). In the modulated case, we observe a fast and long-range spreading of the wavefunction (Fig. 3a), that is responsible for the tremendous growth of the standard deviation (Fig. 3e). As an additionnal signature of the long-range hopping, the standard deviation appears to saturate with a clear finite size effect, that we attribute to the fact that the boundaries are reached very fast. We have also simulated the exact dynamics of the corresponding undriven lattice with no chaos that highlights the clear contrast between the two systems. Indeed, the unmodulated case gives a slow and short-range ballistic spreading of the wavefunction with no finite-size effect (Fig. 3b ande).

Analytical derivation of the hopping law.-In addition to the expected long-range decay ∝ 1/n of the effective coupling term, numerical simulations show apparent erratic fluctuations around this algebraic law (Fig. 4). We discuss here a simple model to explain their origin. For each of the N res resonances in the effective band, we apply a two-level model that involves only three parameters: the slope α = dε ch /dβ of the energy of a chaotic state with β, the coupling intensity W between the chaotic and the regular states and the position β 0 of the crossing in the spectrum (see Fig. 2). Using the linearity of Eq. ( 2) and assuming sharp resonances (∆β 2π/λ), the asymptotic behavior of t eff n is (see Supplemental Material [54])

t eff n ≈ i πn resonances sgn(α)|W |e inβ0λ . ( 3
)
This simple model is in very good agreement with numerical data (see Fig. 4) and shows that the relevant time scale of the tunneling dynamics is heff /|W |. The phase term e inβ0λ , which depends on the position of the resonances in the effective band, gives the observed fluctuations of hopping amplitudes around the algebraic law.

Since the W 's of the N res resonances are associated with tunnel coupling to chaotic states, Random Matrix Theory suggests that they can be described as independent Gaussian variables with a fixed variance w 2 . In the same spirit, as soon as n is large enough the phases nβ 0 λ mod [2π] can be considered random. Using the known results about sums of complex numbers with Gaussian amplitudes and random phases [55], Eq. ( 3) leads to a simple statistical model for the couplings, with |t eff n | ≡ W/n with W a Gaussian random variable of variance N res w 2 . We stress that this implies the distribution of n|t eff n | is universal. Fig. 4b shows the validity of this approach. (∆ being the energy difference with the chaotic state involved). The effective picture is thus legitimated by both (1) the sharpness of the resonances that guarantees that the total part of the system that is delocalized in the chaotic sea is small at any time (the amplitude of oscillations being large only close to the resonances), and (2) on the observation that the slowest Rabi oscillation is from Eq. (3) always faster than the induced tunneling process (h eff /|t eff n | ≥ πh eff /W ). These arguments are corroborated by Fig. 3d: the projection of the system on the chaotic sea displays fast and weak oscillations around a very low value.

The arguments we present emphasize that the longrange property is a direct consequence of the existence of sharp and strong tunneling resonances in the band structure, properties which are fairly generic when the classical dynamics of the driven lattice is mixed.

Conclusion.-In this letter we generalized the original chaos-assisted tunneling mechanism between two wells to spatially periodic lattice systems. We demonstrated that in an intermediate regime of temporal driving, the system dynamics could be mapped to a tight-binding Hamilto-nian with long-range hopping. Strikingly, this new manifestation of chaos at a quantum scale is fairly generic and could be observed in many different experimental situations. In the context of quantum simulation, this result opens the possibility to observe the dynamics of longrange models, and thus to investigate many important phenomena of condensed matter such as glassy physics, many-body localization or quantum multifractality.

FIG. 1 .

 1 FIG. 1. Equivalent representations of chaos-assisted tunneling in a driven lattice. (a) In situ representation: a wavefunction tunnels through different potential wells. (b) Phase space representation: the wavefunction escapes from an island of stable orbits (blue) by regular tunneling, spreads in the chaotic sea (red) and tunnels in another island. (c) Tightbinding representation: the system is composed of N sites, the coupling between i-th and j-th site is proportional to 1/|i-j|.

FIG. 2 .

 2 FIG. 2. Sketch of an avoided crossing between regular |βreg and chaotic |β ch Bloch waves, showing the different quantities that appear in the two-level model: |W | the strength of the coupling between |βreg and |β ch , α the slope of the energy of |β ch and β0 the point of equal mixing. The width of the crossing ∆β is defined as the intersection of the slope of the energies of mixed states at β = β0 with the x-axis. Near β = β0, the eigenstates |β± become a mixture of |βreg and |β ch and form two non-crossing branches. As |β± transforms from one original state to the other, the color code gives the intensity of the mixing through the projection of the eigenstate on the regular state |βreg . Solid black line is the effective regular energy (see text). Husimi representations of selected states are given on top of the stroboscopic phase portrait. Inset: Quasi-energy dispersion relation of the Hamiltonian Eq. (1) (h eff = 0.4, γ = 0.20, = 0.15), black solid line corresponds to the effective regular band and red dashed line to a nearest-neighbor approximation with parameters extracted from the value of the effective regular band at β = 0 and β = π/λ.

24

 24

FIG. 3 .

 3 FIG. 3. Dynamics of a wavepacket, initially located on a single regular site n0. Parameters are γ = 0.2, heff = 0.4 and ε = 0.15 (modulated lattice) or ε = 0 (unmodulated lattice). (a-b) Probability at each site vs time, normalized for visibility at each time by its maximum value over the lattice. Exact dynamics in the modulated (a) and unmodulated (b) lattices and its corresponding effective dynamics (note that the system is symmetric through n -n0 → n0n). (c) Probability at each site for t = 1500T and N = 1079, solid line for exact dynamics and symbols for effective dynamics, red for ε = 0.15 and black for ε = 0. (d) Overlap of the wavefunction with the chaotic sea in the modulated lattice. Same color code as (c), additional blue (N = 539) and green (N = 269) data for modulated lattices of smaller sizes. (e) Standard deviation of the wavefunction (see text), same color code as (d).

FIG. 4 .

 4 FIG. 4. (a) Effective hopping intensity |t eff n | vs distance between sites n for γ = 0.2, ε = 0.15 and heff = 0.4. Blue solid line corresponds to the data extracted from numerical Fourier series of the effective band structure. Red solid line corresponds to Eq. (3) with parameters extracted from the band structure. Black solid line is the typical value of Eq. (3). Inset: same data rescaled to show the small-distance behavior, additional black solid line is the unmodulated case ε = 0. (b) Distribution of fluctuations around the 1/n law for 5 parameter sets: histogram corresponds to cumulative values for 1500 < n < 10000, dots are partial datasets of 500 consecutive values of n, black curve is analytical prediction (see text).

  Discussion.-The theoretical results presented above rely on the effective Hamiltonian picture. It is thus important to assess its validity in our context. The exact tunneling dynamics between two sites can be written(n + m) reg |U F |m reg = 1 N β e iβλn βreg |U F |β reg . In the effective approach β reg |U F |β reg writes exp -iε eff reg (β)t/h eff , thus it does not take into account the Rabi oscillations of each regular Bloch wave |β reg with the chaotic sea |β ch , whose amplitude is given in a two-level approximation by W/ √ W 2 + ∆ 2 and whose period is πh eff / √ W 2 + ∆

  2

This work was supported by Programme Investissements d'Avenir under the program ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT, and research funding Grant No. ANR-17-CE30-0024. We thank Calcul en Midi-Pyrénées (CALMIP) for computational ressources and assistance. We thank O. Gauthé and B. Peaudecerf for useful discussions.