Neutron Stars
Rotational and Magnetic Field Instabilities&
Gravitational Waves

Kostas Kokkotas
Theoretical Astrophysics
Eberhard Karls University of Tübingen

13.05.2016

NAOJ

1

Neutron Stars: Mass vs Radius
Static Models

Rotating Models
supramassive

Demorest etal 2010

𝑀"#$ ≃ (1.1962 + 0.0108)𝑀012
Breu-Rezzolla 2015
20/1/16Frankfurt2

Neutron Stars & “universal relations”
Need for relations between the “observables” and the“fundamentals” of NS physics

20/1/16

Average Density

ρ ~ M / R3

Compactnessz~M Rη = M3 /I

Moment of Inertia

I ∼ MR 2

I ∼ J /ΩQuadrupole Moment

Q ~ R5Ω2

Tidal Love Numbers

λ ~ I 2Q

Frankfurt3

I-Love-Q relations
EOS independent relations were derived by Yagi & Yunes(2013) for non-magnetizedstars in the slow-rotation and small tidal deformation approximations.

… the relations proved to be valid (with appropriate normalizations) even for fast
rotating and magnetized stars
20/1/16Frankfurt

4

NEUTRON STARS &ALTERNATIVE THEORIES OF GRAVITY

13.05.2016

NAOJ

5

STT of gravity - Motivation•

The Scalar Tensor Theory (STT) is one of the most naturalgeneralizations of the Einstein’s Theory of Gravity (ETG)•

Their essence is in one or several scalar fields that are mediators of
the gravitational interaction in addition to the spacetime metric ofclassical ETG•

Scalar fields appear in the reduction of the Kaluza-Klein theories to
4 dimensions, in string theory and in higher dimensional gravity but
STT can be defined completely independently•

STT can be considered as an ETG with variable gravitational
constant•

They fit to the observational data very well•

They are also an essential part of dark energy and dark mattermodels•

The f(R) theories are mathematically equivalent to the STT

13.05.2016

NAOJ

6

II. BASIC EQUATIONSSTT of gravity – Action

The general form of the gravitational action of the scalar-tensor theories in the physicalme is given by [1–4]• Physical (Jordan) frame action:

1
S=
16πG∗

!

4

d x

"

#

µν

$

#

$

− g̃ F (Φ) R̃ − Z (Φ) g̃ ∂µ Φ∂ν Φ − 2U (Φ) + Sm Ψm ; g̃µν .

e, G∗ is the bare gravitational constant, R̃ is the Ricci scalar curvature with respec• Einstein frame action (much simpler):
ce-time metric g̃µν . The dynamics of the scalar field Φ is governed by the functionΦ) and U (Φ). In order for the gravitons to carry positive energy we must have F (Φ) >negativity of the scalar field kinetic energy requires that 2F (Φ) Z (Φ) + 3[dF (Φ)/dΦ
1 2
βϕ on Ψ m amatter fields are collectively denoted by Ψm and theiractionSdependsmd ln ( A(ϕ ))2A(ϕ ) = ek(ϕ ) =
– Coupling
function
ce-time metricg̃µν . Here we
consider the phenomenological
case
when
the matter actiodϕ) =satisfied.
βϕ1involve the scalar field in order for the weak equivalence principlek(
toϕbe
From a mathematical point of view it is convenient to analyze the scalar-tensor theori– We set the potential to zero
V (ϕ ) = 0
pect to the conformally related Einstein frame given by the metric:

13.05.2016gµν =F (Φ) g̃µν .
NAOJ

7STT of gravity – Neutron Stars
Spontaneous Scalarizarion is possible for β<-4.35(Damour+Esposito-Farese 1993)

Properties of the static scalarized neutron stars

The solutions with nontrivial scalar field are energetically more favorable
than their GR counterpart (Harada 1997, Harada 1998, Sotani+Kokkotas 2004).
13.05.2016

NAOJ

8

STT of gravity - Observations|k0 |Observational constraints:

k0<0.004
β> -4.8 (-4.5)
Damour & Esposito-Farese (1996, 98)
Will (2006), Freire etal (2012)

Freire at al (2012)13.05.2016

NAOJ

9

STT of gravity – Fast Rotating Stars•

The effect of scalarization is much stronger for fast rotation.

•

Scalarized solutions exist for a much larger range of parameters than in thestatic case
Doneva, Yazadjiev, Stergioulas, Kokkotas 201313.05.2016

NAOJ

10

STT of gravity
Angular Momentum & Moment of Inertia

Models with constant central energy density

Sequences of models rotating at the Kepler
limit

Not surprizing that both angular momentum and moment of inertial could differ
twice for scalarized solutions
13.05.2016

NAOJ

11

NSs in f(R)-gravity: Static Models
f (R) = R + aR••

2

The differences between the R2 and GR are comparable with the uncertainties in the
nuclear matter equations of state.
The current observations of the NS masses and radii alone can not put constraints on
the value of the parameters a, unless the EoS is better constrained in the future.
Yazadjiev, Doneva, Kokkotas, Staykov (2014)13.05.2016

NAOJ

12

NSs in f(R)-gravity: Fast Rotationf (R) = R + aR 2

Mass of radius diagrams for two realistic EOS

Difficult to set constraints on the f (R) theories using measurement of the neutronstar M and R alone, until the EOS can be determined with smaller uncertainty.
Yazadjiev, Doneva, Kokkotas, (2015)13.05.2016

NAOJ

13

NSs in f(R)-gravity: Fast Rotationf (R) = R + aR 2

Yazadjiev, Doneva, Kokkotas (2015)

ü The differences in the neutron star moment of inertia on the other hand can be much
more dramatic.
ü Large deviations can be potentially measured by the forthcoming observations of the NS
moment of inertia [Lattimer-Schutz 2005, Kramer-Wex 2009] that can lead to a direct testof the R2 gravity.
13.05.2016

NAOJ

14

NSs in f(R)-gravity: I-Q relations / Fast Rotationf (R) = R + aR 2
I ≡ I /M3
Q ≡ Q / (M 3 χ 2 )χ ≡ J /M2
ν : rot. frequency (Hz)≡Doneva, Yazadjiev, Kokkotas (2015)• The results show that the I-Q relation remain nearly EoS independent for fixedvalues of the normalized rotational parameter• The differences with the pure Einstein’s theory can be large reaching above 20%
for lower masses and slow rotation.
13.05.2016

NAOJ

15

Oscillations & Instabilities
The most promising strategy for constraining the physics of neutron stars involves
observing their “ringing” (oscillation modes)
Ø f-mode : scales with average density
Ø p-modes: probes the sound speed through out the star
Ø g-modes : sensitive to thermal/composition gradients
Ø w-modes: oscillations of spacetime itself.
Ø s-modes: Shear waves in the crustØ Alfvèn modes: due to magnetic fieldØ i-modes: inertial modes associated with rotation (r-mode)
Typically SMALL AMPLITUDE oscillations –> weak emission of GWs
UNLESS
they become unstable due to rotation (r-mode & f-mode)𝑙 = 2, 𝑚 = 2

13.05.2016

𝑙 = 3, 𝑚 = 3

𝑙 = 4, 𝑚 = 4

NAOJ

16

MAGNETARS:A PROMISING CASE FOR ASTEROSEISMOLOGY

13.05.2016

NAOJ

17

Magnetars•

Young, slowly spinning (P~10s) systems (20+)•

Exhibit regular γ-ray flares– Believed to be powered by magnetic field– Either trigger or are preceded by starquakes– Some linked to glitches or anti-glitches•

Three giant flares observed with peakluminosities ~1047 erg/s
–
–
–
–

•

March 5, 1979 :SGR 0526-66
August 27, 1998 :SGR 1900+14
December 27, 2004:
SGR 1806-20
Recently few medium ones

Giant flares– QPOs – 10’s -100’s of Hz– Magnetic field reconstruction

13.05.2016

NAOJ

18

Magnetars:
Quasi-Periodic Oscillations (QPOs)
ü Giant flares in SGRs•

A decaying tail for several hundred seconds follows the flare.

ü QPOs in decaying tail (Israel et al. 2005; Watts & Strohmayer 2005, 2006)•

SGR 1900+14 : 28, 54, 84, and 155 Hz•

SGR 1806-20 : 18, 26, 29, 92.5, 150, 626.5, 720, 976, 1837, 2384 Hz•

SGR 1806-20 : Additional frequencies 22, 16, 60, 116 Hz, also 720 & 2384 Hz;(Hambaryan, Neuhaeuser, Kokkotas 2011)

QPOs

13.05.2016

NAOJ

19

Alfven Continuum and/or Discrete oscillations
12

Only Crust Oscillations

Sotani, Kokkotaa, Stergioulas 2007,2008
Samuelsson, Andersson 2007
Sotani, Colaiuda, Kokkotas 2008
Steiner, Watts 2009…
Sotani etal 2012-16

10

Polytropic+NV14 B=1016 l=8 210 Hz

8
y

•
•
•
•
•
•

6

Without Crust••••

Levin 2007
Sotani, Kokkotas, Stergioulas 2008
Colaiuda, Beyer, Kokkotas 2009
Cerda-Duran, Stergioulas, Font 200942

Fluid + Crust
•
•
•
•
•

Van Hoven, Levin 2011, 2012
Cerda-Duran, Stergioulas, Font 2011
Colaiuda, Kokkotas 2011
Gabler etal 2012
Gabler etal 2013 …

Superfluidity• Passamonti, Lander 2012
• Sotani etal 2013
• Gabler etal 2013

Mixed axial-polar• Colaiuda, Kokkotas 2012
• Lee, Yoshida 2015

Non-axisymmetric• Sotani, Kokkotas 201213.05.2016

0

2

4

6

8

10

12

x

ü The combination
Ø of poloidal + toroidal magnetic fields+ crustØ leads to coupling between polar andaxial modes leading toØ PURE discrete spectrum
ü The main results of the magnetar
seismology remain unchanged !
NAOJ(Colaiuda-Kokkotas 2012)20

Magnetars: SGR 1806-20
2005-612Israel etal
Watts &
Strohmayer

Polytropic+NV14 B=1016 l=8 210 Hz

10

y

8

EoS :
Mass:

6
4

2
crustal modes

4
APR 1418

B = 4 × 10 15

Gauss

53

66

16

16

2013-15

18

26

29.5

30
39

22

37

44

discrete Alfvén modes

22 26

Hambaryan,
Neuhaüser, KK

26

x
Crust
: 0.099 R

2011

ColaiudaKK

22

Radius: 11.57 km
15 G
B-field: 2x10
6
8
10
12

2

0

APR (NV)
M=1.4M¤

2010-12

53

82

56
60

56
59 &61

66
74
16 18

30

39

44

56 60
Frequency(Hz)

74

80

80

92

82

Colaiuda, Kokkotas (2011-12)
13.05.2016

NAOJ

92

92

92

21

Magnetars: Open Questions
Do we understand how the QPOs are excited?
The answer is NO! (at least partially)
Great progress in the last 7-8 yearsü BUT mainly AXISYMMETRIC oscillations used to explain the observed
QPOs
ü For NON-Axisymmetric oscillations both poloidal & toroidal B-fields
are unstable!

A new event of the type of SGR 1806-20, might be catalytic for understanding:• The mechanism that triggers the hyperflares• The QPOs in the decaying tail• The EOS, the Mass, Radius, B-field of magnetars
13.05.2016

NAOJ

22

1973M

The Tayler Instability
Toroidal Fields
Bernstein etal (1958) showed that the stability of an ideally conducting
system depends on the sign of the change of potential energy of the system
for an arbitrary perturbation ξ(x,t):

 
 

1 % 2   
δW = ∫ &Q − j ⋅ Q × ξ + γ P(divξ ) + ξ ⋅ gradPdivξ + ξ ⋅ gradΦdivρξ '(dτ2 Q = curl ξ × B()

For a toroidal field the m=0 and the m=1 will look like:
Rotation or/and a strong poloidal field can potentiallywork against this instability

13.05.2016

NAOJ

23

Magnetic Field Stability: Semianalytic Approach
Herbrick+KK 2015

δ W = ∫ ⎡⎣δ WB + δ W fluid + δ Wgrav ⎤⎦ dV

In agreement with Akgün etal (2013)V

ü Pure poloidal
ü Pure toroidal
ü Mixed field + stratification

The pure toroidal field would beunstable but can be stabilized by thecontribution from the poloidal field thatis positive in this area.

ü Constraining realistic magnetic field structures by parametrisation
ü Adding more realistic features e.g. crust, superconductivity,…
ü Constructing arbitrary displacement fields
ü Eventually providing a GR criterion of stability
13.05.2016

NAOJ

24

Simulation of Magnetic Field Instability
Lasky, Zink, Kokkotas, Glampedakis ApJL (2011)13.05.2016

NAOJ

25

Simulation of Magnetic Field Instability
Lasky, Zink, Kokkotas, Glampedakis ApJL (2011)13.05.2016

NAOJ

26

Gravitational Waves from Magnetars
Zink, Lasky, KK (2012)13.05.2016

NAOJ

27

Gravitational Wave Asteroseismology
Oscillation patterns can reveal theinternal structure of neutron stars :ü mass,ü radius,
ü EoS,
ü rotation,
ü B-field,
ü crust,…
4.0

Andersson,Kokkotas 1996,1998,2001

150
A
B
C
D
E
F
G
I
L
G240
G300
WFF
Q_1
Q_2
Q_3

140

3.6

130

3.2

2.8

2.4

2.0

1.6

1.2

0.03

0.04

0.05

0.06

3 1/2

13.05.2016(M/R )

0.07

0.08

120

Rωw -modeωf-mode

A
B
C
D
E
F
G
I
L
G240
G300
WFF
Q_1
Q_2
Q_3

110

+

100
90
80

You getM/R & M/R3

70
60
50

0.15

0.18

0.21

0.24

M/R

0.27

NAOJ

0.30

Lattimer+Prakash 200728

f-modes: Asteroseismology
We can produce empirical relation relating the parameters of therotating neutron stars to the observed frequencies.
Gaertig-Kokkotas 2008, 2010, 2011

Frequency
Cowling Approximation

Damping/Growth time

13.05.2016

NAOJ

29

Asteroseismology: Realistic EoS
Doneva, Gaertig, KK, Krüger (2013)

Nearly “universal” fitting formulae for :• the frequencies• the damping times• Independent of GR or Cowling⎛ ωc ⎞⎛ Ω ⎞≈ f⎜⎜⎝ ω ⎟⎠⎝ Ω K ⎟⎠
0 ℓ=2,3,4
1.2⎛ τ0 ⎞⎜⎝ ⎟⎠τm)ranch (l=unstable b

1.0
stab

le b
r

ω/ω0

0.8
0.6
0.4
0.0

anc
h

(l=m

1/2

⎛ω ⎞
≈ f⎜ i⎟⎝ ω0 ⎠

)

l=m=2
l=m=3
l=m=4
l=-m=2, 3, 4
l=2 Full GR, C and S models
l=3 Full GR, C and S models

0.2

0.4

0.6

0.8

1.0

Ω/Ωk

Oscillation frequencies
13.05.2016

Damping/Growth Times
NAOJ

30

Asteroseismology
Stable Branch

Unstable Branch

Unstable Branch

Doneva, Gaertig, KK, Krüger (2013)13.05.2016

NAOJ

31

Asteroseismology: alternative scalings
M σ iunst = [ (0.56 − 0.94ℓ) + (0.08 − 0.19ℓ)M Ω + 1.2(ℓ + 1)η ]
14

l=m=2
l=m=3
l=m=4

- b(M

)

12
10
8
6
4
1.2

The l = 2 f-mode oscillation frequencies
as functions of the parameter η

13.05.2016

NAOJ

1.4

1.6

1.8

2.0

2.2

η = M3 /I
Doneva-Kokkotas 201532

Asteroseismology: alternative scalings⎛ M ⎞The normalized damping time η ⎜⎝ τη 2 ⎟⎠(1/2ℓ)whereη = M3 /I

as a function of the normalized oscillation frequency Mσ for l = m = 2 & l = m = 4 fmodes.
Doneva-KK 201513.05.2016

NAOJ

33

Asteroseismology:
Alternative Theories of Gravity•

The maximum deviation between the f-mode frequencies in GR and R 2 gravity
is up to 10% and depends on the value of the R 2 gravity parameter a .

•Alternative normalizations show nicer relations

η = M3 /I
13.05.2016

NAOJ

34

The CFS instability
Chandrasekhar 1970: Gravitational waves lead to a secular instability
Friedman & Schutz 1978: The instability is generic, modes with sufficiently large m are
unstable.

A neutral mode of oscillation signals the onset of CFS instability.
ü

Radiation drives a mode unstable if themode pattern moves backwards
according to an observer on the star
(Jrot<0), but forwards according tosomeone far away (Jrot>0).

üThey radiate positive angularmomentum, thus in the rotating frame
the angular momentum of the modeincreases leading to an increase inmode’s amplitude.

Gaertig+Kokkotas 2008ωωin= − rot + Ωmm
LIGO/Virgo/GEO/KAGRA/ET band
13.05.2016

NAOJ

35

The Excitation of Secular Instabilities
SATURATION AMPLITUDE REACHEDOSCILLATION AMPLITUDE

Mode is damped via shocks ormode coupling
UNSTABLE GROWING PHASE
Duration a few sec – hours

Rotation
Temperature
EoS,…

UNKNOWNS(?)üCritical rotationüSaturation AmplitudeüDurationüDifferential RotationüInstability WindowüInfluence of B-fields

STABLE PHASE

TIME
13.05.2016

NAOJ

36

The CFS Instability Window
1.00

Kepler limit

f-modeMutual friction

ROTATION Ω/ΩK

Lindblom&Mendel‘95

0.75
Shear
Viscosity
Magnetic field
Hyperon viscosity
Crust
Quark matterr-modeBulk
Viscosity

0.05

107 K

108 K

109 K

1010 KTEMPERATURE
13.05.2016

NAOJ

37

Instability Windowüü

For the first time we have the window of f-mode instability in GR
Newtonian: (l=m=4) Ipser-Lindblom (1991)>30 min

Mutual friction

N=0.6613.05.2016

NAOJ

Gaertig, Glampedakis, Kokkotas, Zink (2011)38

Saturation of the Instability
Parametric Resonance13.05.2016

NAOJ

Pnigouras, Kokkotas (2015)39

Saturation of the Instability
Parametric Resonance13.05.2016

NAOJ

Pnigouras, Kokkotas (2015)40

Saturation of the Instability
Parametric Resonance13.05.2016

NAOJ

Pnigouras, Kokkotas (2015)41

Evolution of a nascent (unstable) NS

Mutual Friction plays NO ROLE for the f-mode instability
Procedure as described in Owen etal 1998 & Anderson, Jones, KK 2002
Passamonti-Gaertig-KK-Doneva (2013)13.05.2016

NAOJ

42

WFF2 Mb= 1.8
f

S/N : aLIGO

3.5=10-6, l=m=3

d = 10 Mpc
aLIGO
ET

3.0

40
35
30

2.5

25

2.0

20

1.5

15

1.0

10

0.5

5

0.0

1010

1011

B

1012

1013

2.0

WFF2 Mb= 1.8
f0

d = 10 Mpc
aLIGO
ET

1.5

20

15

1.0

10

0.5

5

0.0

1014

=10-7, l=m=3

1010

1011

B

1012

1013

S/N : Einstein Telescope

4.0

45

S/N : aLIGO

4.5

S/N : Einstein Telescope

Evolution of a nascent (unstable) NS

0

1014

The instability can be potentially observed by events in Virgo cluster
BUT
•
•
•

Event rate is unknown
Competition with r-mode and magnetic field slow-downSaturation amplitude is varying during the procces
Passamonti-Gaertig-Kokkotas-Doneva (2013)13.05.2016

NAOJ

43

A GRAVITATIONAL WAVE AFTERGLOWIN BINARY NEUTRON STAR MERGERS

13.05.2016

NAOJ

44

Binary Neutron Star Mergers
the standard scenario
I.

After the merging the final body mostprobably will be a supramassive NS (2.53 M¤)

h+,× [10-22]

2 (c) H15
1
0

-1
-2

The body will be differentially rotatingh+,× [10-22]

III. The “averaged” magnetic field will
amplified due to magnetic fieldinstabilities (up to 3-4 orders ofmagnitude)

1
0

-1
-2

-5

This phase will last only a few tenths ofmsecs and can potentially provide
information for the Equation of State(EOS)

h+,× [10-22]

2 (a) H135

IV. The strong magnetic field and theemission of GWs will drain rotationalenergyV.

2 (d) S15

NAOJ

5

10

15

5

10

15

tret - tmerge [ms]10

-2
2 (b) S135
1
0

-1
-2

-5

13.05.2016

0-1h+,× [10-22]II.

0tret - tmerge [ms]Kiuchi, Sekiguchi, Kyutoku, Shibata2012
45

Post-Merger Scenario
Three different outcomes of the merger of a BNS merger
Rowlinson 2013

ü The outcome is dependent upon the mass (M) of the central object formed and the
maximum possible mass of a neutron star (Mmax).
ü On the right are sketches of the expected light-curves if a stable (top) or an
unstable magnetar (bottom) is formed.
13.05.2016

NAOJ

46

Short γ-ray light curves§

The favored progenitor model for SGRBs is the merger of two NSs that triggers an explosion
with a burst of collimated γ-rays.

§

Following the initial prompt emission, some SGRBs exhibit a plateau phase in their X-ray light
curves that indicates additional energy injection from a central engine, believed to be a rapidly
rotating, highly magnetized neutron star.

§

The collapse of this “protomagnetar” to a black hole is likely to be responsible for a steep decay
in X-ray flux observed at the end of the plateau.

Rowlinson, O’Brien, Metger,
Tanvir, Levan 201313.05.2016

NAOJ

47

Post-Merger NS: secular instability
Doneva-KK-Pnigouras 2015
The post-merger object is still stable and rotates at nearly Kepler periods < 1ms
The evolution into the instability window
NS*evolu3onary*path*α (Mc2 )
10-2
10-4
10-6
Black&Hole*Limit*

10-8
10-10
10-12
10-14

The detailed evolution depends:

13.05.2016a) Strength of the magnetic field (averaged mayreach 1015-16 G !)b) Equation of state of the post-merger neutron starc) Fine details of the non-linear dynamics (threemode coupling, shock waves, wave breaking)
NAOJ
48

S/N : aLIGO

f=10-5, l=m=2

d = 50 Mpc
aLIGO
ET

3.0
2.5 APR Mb= 3.2
2.0

30
25
20

WFF2 Mb= 3.0

1.5

15

1.0

10

0.5
0.0

1012

1013

B

1014

1015

4.0

f

3.0

=10-7, l=m=2

d = 50 Mpc
aLIGO
ET

2.5 APR Mb= 3.2
2.0
1.5

13.05.2016

B

2.0
1.5

25
20

WFF2 Mb= 3.0

15

0.0

35
30
25

5

1015

2.5 APR Mb= 3.2

30

0

0.5

1014

d = 50 Mpc
aLIGO
ET

3.0

35

5

15

1013

=10-6, l=m=2

0.5

10

1012

f

5

1.0

0.0

WFF2 Mb= 2.9

10

20
WFF2 Mb= 3.0

3.5

40

1.0

40

3.5 WFF2 Mb= 2.9

S/N : aLIGO

35

S/N : aLIGO

WFF2 Mb= 2.9

3.5

4.0

S/N : Einstein Telescope

40

1012

1013

B

1014

1015

0

Competition between the B-field andthe secular instability

S/N : Einstein Telescope

4.0

S/N : Einstein Telescope

Post-Merger NS: secular instability

GW frequencies:
WW2a: 920-1000 Hz
APR: 370–810 Hz
WFF2b: 600–780 Hz
Doneva-KK-Pnigouras 2015

0

NAOJ

49

Conclusions
ü The influence of alternative/extended theories of gravity on
NS parameters is much more pronounced for fast rotation.
ü Difficult to set constraints on theories using measurementof the neutron star M and R alone, until the EOS can be
determined with smaller uncertainty.

13.05.2016

NAOJ

50

Conclusions
ü The influence of alternative/extended theories of gravity on
NS parameters is much more pronounced for fast rotation.
ü Difficult to set constraints on theories using measurementof the neutron star M and R alone, until the EOS can be
determined with smaller uncertainty.
ü Asteroseismology for fast rotating stars is possible
ü Asteroseismology for magnetars is promising (!)

13.05.2016

NAOJ

51

Conclusions
ü The influence of alternative/extended theories of gravity on
NS parameters is much more pronounced for fast rotation.
ü Difficult to set constraints on theories using measurementof the neutron star M and R alone, until the EOS can be
determined with smaller uncertainty.
ü Asteroseismology for fast rotating stars is possible
ü Asteroseismology for magnetars is promisingü f-mode instability can be potentially a good source for GWs
for supramassive NS
ü The efficincy depends on the saturation amplitude andstrength of B-field.

13.05.2016

NAOJ

52